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ABSTRACT 

High stagnation pressures and enthalpies are required for the testing of aerospace 
vehicles such as aerospace planes, aeroassist vehicles and reentry vehicles. Among the 
most useful ground test facilities for performing such tests are shock tunnels. With a 
given driver gas condition, the enthalpy and pressure in the driven tube nozzle reservoir 
condition can be varied by changing the driven tube geometry and initial gas fill pressure. 
Reducing the driven tube diameter yields only very modest increases in reservoir pressure 
and enthalpy. Reducing the driven tube initial gas fill pressure can increase the reservoir 
enthalpy significantly, but at the cost of reduced reservoir pressure and useful test time. A 
new technique, the insertion of a converging section in the driven tube is found to produce 
substantial increases in both reservoir pressure and enthalpy. Using a one- dimensional 
inviscid full kinetics code, a number of different locations and shapes for the converging 
driven tube section were studied and the best cases found. For these best cases, for driven 
tube diameter reductions of factors of 2 and 3, the reservoir pressure can be increased by 
factors of 2.1 and 3.2, respectively and simultaneously the enthalpy can be increased bj 
factors of 1.5 and 2.1, respectively. 


I. INTRODUCTION - THE REFLECTED SHOCK SHOCK TUNNEL 

The basic setup of a reflected shock shock tunnel is shown in Fig. 1. It comprises a 

driver tube, driven tube, nozzle and test section. Typically, the driver and driven tubes are 
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separated by a heavy diaphragm and the driven tube and nozzle are separated by a very 
light diaphragm. The driver is initially filled with a high pressure, high sound speed gas, 
typically helium or hydrogen, which may be heated to further increase its sound speed and 
enthalpy. The driven tube is initially filled with the test gas, frequently air or nitrogen, at 
a much lower pressure. Upon rupture of the main diaphragm a shock wave travels down 
the driven tube and an expansion wave system moves into the driver gas. The shock wave 
reflects from the driven tube end wall and breaks the light secondary diaphragm. The 
doubly shocked gas forms a nozzle plenum reservoir condition. Gas from this reservoir 
rushes through the nozzle and flow is established there and in the test section. The useful 
test time is typically hundreds of microseconds to 5 or 10 milliseconds. 

We now examine the shock tunnel operation in more detail using the “X-T diagram 
shown in Fig. 2. In Fig. 2, X is the distance along the driver and driven tubes and T is 
the time after the rupture of the diaphragm. Figure 2 is schematic only and is not to scale. 
In Fig. 2, shock waves are shown as heavy lines, undefined (compression or expansion) 
waves as single light lines, expansion wave systems as fan-like systems of lines and the 
interface between the driver and driven gas as a dashed line. The numbers within the 
zones identify regions of various thermodynamic states. The initial states of the driver and 
driven gas are 4 and 1, respectively. Upon rupture of the diaphragm, shock Si moves down 
the driven tube, compressing the gas from state 1 to state 2. SI reflects at the driven tube 
end wall as shock S2, compressing the gas from state 2 to state 5. The driver gas expands 
through two wave systems, El and E2, from state 4 to state 3. The gas in state 3 drives 
the initial shock SI through the driven gas. In the classic constant area shock tube the 
expansion systems El and E2 become one continuous expansion system and the zone 3 
does not exist. With an area contraction between the driver and driven tubes, the part of 
the unsteady expansion wave system between lines P and Q is replaced by a quasi-steady 
expansion through the area contraction between states 3’ and 3” and an extended region 
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at state 3’ appears. 

The reflected shock S2 passes through the interface at M. If the acoustic impedances of 
regions 2 and 3 are properly matched, there will be no reflected wave W, all the additional 
waves near the interface above M will likewise vanish and regions 5, 6, 7, 8, etc. become one 
region at state 5. This is called “tailored interface” operation. Tailoring can be achieved 
by varying the initial pressure in region 1. The plenum reservoir for the nozzle during the 
test time is then at state 5. If one is not operating at the tailored interface condition, wave 
W and the additional reflected waves above W will exist and conditions in regions 5, 6, 7, 
8', etc., will be different. However, in practice, after 2 or 3 reflections from the interface, 
the waves become quite weak and regions 7, 8, and beyond are essentially at the same 
condition. This condition will be essentially at the pressure which would be generated by 
a single shock propagating into region 3, bringing that gas to “rest . This is referred to as 
“equilibrium interface” operation. The plenum reservoir for the nozzle during the test time 
is then at the state of regions 7, 8, etc. (which are very closely at the same thermodynamic 
state). 

The available test time of shock tunnels is limited by several factors. Arrival of the 
head of the reflection of the driver expansion wave system (HE) at the nozzle entrance 
at time TE initiates a rapid drop in nozzle reservoir pressure and will end the useful test 
time. Also, passage of the driven gas-driver gas interface I through the nozzle entrance at 
time TI will end the useful test time. In equilibrium interface operation, the time for the 
nozzle reservoir condition to settle out, say, between TSH and TS, is not useful test time. 
Although not related directly to Fig. 2, time is required for the establishment of steady 
flow in the nozzle and test section and this time cannot be use for testing. 

The interface motion shown in Fig. 2 assumes a perfect one-dimensional interface. As 
is well known, this is very far from being the case in real shock tubes and shock tunnels. 
A number of effects are known to cause substantial spreading out of the interface. These 
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include interface distortion upon diaphragm opening, boundary layer effects, 1 3 bifurcation 
of the reflected shock (S2) near the wall, 4 instability of the interface 2 - 3 ’ 5-7 and combustion 
at the interface. 3 We do not go into detailed discussion of these effects here, but note the 
following. For incident shock Mach numbers of 8-10, Fig. 15 of Ref. 3 gives a correlation 
of experimental data for driver gas free test times. These are only ~25% of the perfect 
one-dimensional interface values. Also, for some conditions, e.g., some data points of Fig. 
14 of Ref. 3, the driver gas free test time can drop to nearly zero. 

There is a continuing need for ground facilities with higher stagnation pressures and 
enthalpies to allow closer simulation in the testing of current and proposed aerospace 
vehicles. The vehicles include the national aerospace plane, aeroassist space vehicles and 
reentry vehicles for earth and other planets. The shock tunnel is one of the types of facilities 
which can provide these simulations. A number of techniques can be used to generate the 
high enthalpy driver gas needed for such shock tunnels. The driver gas can be characterized 
by its specific heat ratio (7), sound speed (c 0 ), enthalpy (h) and escape speed (u e ). The 
last speed is the maximum velocity the driver gas would reach if it expanded into vacuum 
and gives a rough estimate of the maximum initial shock velocity which could be achieved 
if the inital driver gas to driven gas pressure ratio were very large. For an ideal driver 
gas u e can be shown 8 to equal 2c 0 /(7-l). Low molecular weight driver gases, hydrogen 
and helium, are the obvious choices to yield high driver gas sound speeds and enthalpies. 
These driver gas sound speeds can be further increased by heating. The driver gases 
have been heated by electrical resistance heaters, 9 combustion of hydrogen and oxygen, 10 
piston compression, 11 electric arcs, 12,13 and high explosives. 14 Table 1 shows representative 
performances for the driver gases for the first three of these techniques, as well as that 
of room temperature helium and hydrogen. The numbers shown in Table 1 are based on 
ideal gas calculations and are therefore estimates, not exact numbers. The heated driver 
gases in the last three rows can produce the most useful simulations of aerospace vehicles. 
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Table 1 Performances of various driver gases 


Reference 

Gas 

T 

(K) 

7 

c 0 

(km/sec) 

U e 

(km/sec) 

h 

(J/gm) 

— 

He 

300 

1.67 

1.02 

3.04 

1550 

— 

h 2 

300 

1.40 

1.32 

6.61 

4370 

9 

h 2 

700 

1.40 

2.02 

10.10 

10200 

10 

He/H 2 0 

— 

1.50 

2.13 

8.54 

9090 

11 

He 

4660 

1.67 

3.99 

11.90 

24100 


The electric arc shock tube of Refs. 12 and 13 can produce very high shock velocities, 
up to 50 km/sec, but has the disadvantages of a driver pressure limitation of 500 atm, very 
short test times (2-20 ^sec) and large capacitive energy storage requirements (1.24 MJ). 
Because of the pressure and test time limitations, it is not suitable for many of the tests 
needed for aerospace vehicles. The high explosive driven shock tube of Ref. 14 destroys a 
large amount of hardware with each test run and cannot be considered a facility useful for 
an extensive research or test program. Hence, we will focus our attention on the electrical 
resistance heated driver, the combustion heated driver and the piston compression heated 
driver of Refs. 9-11. 

II. VARIATION OF DRIVEN TUBE DIAMETER AND INITIAL FILL PRESSURE 

Let us consider a given driver performance, i.e., with sound speed and enthalpy given, 
say, by one of the last three rows of Table 1. There will, of course, be a pressure limitation 
also for the driver. The nozzle plenum reservoir pressure and enthalpy can be changed 
by changing the driven tube area and initial fill pressure. We will discuss these variations 
with reference to the NASA Ames 16 Inch Shock Tunnel. The driver tube is 21 m long, 43 
cm inside diameter and the driven tube is 26 m long, 30 cm inside diameter. The nozzle 
is 5.9 m long, the exit diameter is 99 cm and the area ratio can be varied from 95 to 271. 
The pressure ratings of the driver and driven tubes are 680 atm. The combustion driver 
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operates with a mixture of 2H 2 + 1.170 2 +9.36He ignited by 4 wires strung the length of 
the tube and heated by a capacitor bank discharge. Further details can be found in Ref. 
10 . 

The performance of the Ames Shock Tunnel was analyzed using a simple zoned, invis- 
cid, ideal gas computational method described in Appendix A. This method does not give 
exact predictions of tunnel performance but does show the trends of the effects of varying 
the driven tube diameter and initial fill pressure. There are only a few percentage points 
differences between the enthalpies and pressures reached in the nozzle reservoir conditions, 
depending on whether the nozzle is plugged or not. Hence, because the calculations are 
simpler, we consider below only cases with the nozzle plugged. The calculations were done 
for 2H 2 + 0 2 + 9He as the drive gas (which is very close to that actually used) and air as 
the driven gas. Figure 3 shows the ratio of tailored interface driver tube reservoir pressure 
to driver after burn pressure (Pres/Pdr) versus the ratio of driven tube area to driver tube 
area (A^n/A^r). The point is an experimental value from the Ames Shock Tunnel. The 
maximum pressure gain from Adn/Adr going from 1 to 0 is about a factor of two. For 
a reasonable value of A^n/A^r, say, 0.3, the gain over a constant area tube is factor of 
~1.5. There will be a small attendent gain in reservoir enthalpy, of the order of 10%, also. 
Making the driven tube smaller has the disadvantages that the amount of gas in the nozzle 
reservoir condition is reduced and boundary layer effects in the driver tube become more 
severe. Clearly, the increases in nozzle reservoir pressure and enthalpy obtainable in this 
way are very limited. 

A second approach is to operate the tunnel in the equilibrium interface mode. Figure 
4 shows calculated results for equilibrium interface operation of the Ames Shock Tunnel 
for the actual value of A dn /A dr , 0.498. Again, the calculations were made for the plugged 
nozzle condition. The abscissa is the ratio of the initial driven tube fill pressure to that at 
the tailored interface condition. The R w and R p curves give the ratios of the driven tube 
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reservoir enthalpies and pressures achieved to those at the tailored interface condition. 
The Rr curve gives the ratio of the compression ratio of the driven gas between the initial 
fill condition and the reservoir condition to that which occurs at the tailored interface 
condition. The 1,1 coordinate point in Fig. 4 is the tailored interface condition. First, 
we see that very little increase in driven tube reservoir pressure can be obtained using 
equilibrium interface operation, perhaps 10%, and this is achieved at at a cost of significant 
drops in enthalpy. 

By decreasing the driven tube initial fill pressure below the tailored interface value, 
significant increases in enthalpy can be achieved. However, these are achieved at a cost 
of decreased reservoir pressure and test time. The decreased reservoir pressure can be 
obtained directly from Fig. 4. Also, however, the driven tube gas is more highly compressed 
and at a higher sound speed than for the tailored interface case. Both these factors will tend 
to make the slug of compressed driven tube gas pass through the nozzle more quickly than 
for the tailored interface case. Further, for the equilibrium interface case, one must wait 
for the reservoir condition to settle out (i.e., wait between TSH and TS in Fig. 2) before 
the test time can begin. This waiting is not necessary for tailored interface operation. For 
example, from Fig. 4, if the initial driven tube fill pressure is reduced to 0.22 times the 
value for tailored interface operation, the reservoir enthalpy will be doubled, the reservoir 
pressure will be reduced 21% and the gas compression ratio will be increased by a factor of 
1.-79. This would cause the time for test gas slug to pass through the nozzle to be reduced 
to about 40% of the value for tailored interface operation. A further reduction in test time 
would occur due to the required wait for settling out of the reservoir pressure. Summing 
up, equilibrium operation can produce only very minor increases in driven tube reservoir 
pressure. It can produce significant increases in driven tube reservoir enthalpy, but at the 
cost of significant reductions in reservoir pressure and test time. 

III. NEW CONVERGING DRIVEN TUBE TECHNIQUE 
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A. Technique 

As was pointed out in the previous section, simply decreasing the size of the driven 
tube with respect to the driver tube produces only very limited increases in reservoir 
pressure pressure and enthalpy. The situation is entirely different if a portion of the driven 
tube is made converging. We will demonstrate that large increases in reservoir pressure 
and enthalpy can be obtained if a converging driven tube is used. A number of geometries 
studied herein are shown in Fig. 5. The diameters of the driver and driven tube are in the 
correct proportion in the figure but the length- to-diameter of the driven tube is not. As 
can be seen from Fig. 5, the converging section of the driven tube can be at the upstream 
end of the tube (cases d,e,f,g,h,j), the downstream end of the tube (cases b,i) or the entire 
driven tube can be converging (case c). Diameter ratios of 2 (most cases) and 3 (cases i,j) 
were studied. The majority of the converging sections had a linear diameter taper, but, in 
addition, both trumpet (cases f,h) and bell shapes (case g) were studied. The equations 
of the driven tube radius ( a taper equations”) as functions of distance along the tube are 


given below. 


Linear taper, case d 
Bell, case g 
Trumpet, case f 
Trumpet, case h 


D x 

— = 1-0.5 — 
Dn Ax 


— = Jl -0.75-^- 
D 0 V Ax 


— = 1 -0.5 A( a 
D„ V Ax 


— = 1 -0.647,/ — 


D 0 


Ax 


(1) 

( 2 ) 

( 3 ) 

( 4 ) 


where D 0 is the diameter of the driven tube inlet, D is the diameter of the driven tube 
at distance x from the inlet and Ax is the length of the contracting section of the driven 
tube. An abbreviated numerical code describing each case is also given in Fig. 5. The first 


digit in each code is the total diameter contraction ratio of the driven tube relative to the 
reference, constant area driven tube. If there is a single number after the comma, it denotes 
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either a constant area driven tube ( a 0” after the comma) or that the contraction occupies 
the entire driven tube (“1” after the comma). Two numbers after the comma denote the 
location in the driven tube where the contraction takes place, as follows. “2,20’ denotes 
that the contraction takes place in the upstream half of the tube, while 2,02 denotes 
that the contraction takes place in the downstream half of the tube. “2,40” denotes that 
the contraction takes place in the upstream quarter of the tube. If there are only numbers 
in the designation, the taper is linear in diameter. Letters at the end of the designation 
denote other curved taper shapes. These include trumpet shapes ( tl , Fig. 5f, t2 , Fig. 
5h; “t3”, to be discussed at a later point), bell shapes (“b”, Fig. 5g) or hyperbolic tangent 
shapes (“h”, not shown in Fig. 5; will be discussed at a later point). 

B. Run Conditions and Computational Method 

The conditions in the driver for all cases studied were calculated as follows. Before 
combustion, the driver was assumed to be filled with a 2 H 2 + 1.1702 + 9.36He mixture 
at 8.622 x 10 7 d/cm 2 and 292 Iv. This gas was then allowed to burn adiabatically to 
an equilibrium mixture and the enthalpy noted. This calculation was repeated with an 
enthalpy loss of 2.5% added in to bring the calculated sound speed into good agreement 
with a representative experimental value of 2.13 km/sec. The resulting burned driver gas 
was at a pressure of 7.075 x 10 8 d/cm 2 , a temperature of 2631 K and had a molecular weight 
of 6.839. The driven gas for most runs was air at 8.67 x 10 5 d/cm 2 and 295 K. The nozzle 
from the secondary diaphragm to the throat is initally filled withe air at 1.014 x 10 2 d/cm 2 
and 295 K. These conditions are nearly identical to those of a number of experimental runs 
made recently (1991) in the Ames 16 In. Hypersonic Shock Tunnel. The only difference is 
that the theoretical cases studied herein are at the maximum rated pressure of the Ames 
facility which is about 25% higher than the pressure of any experimental run made to date. 

The calculations were done using a one dimensional, inviscid CFD code. The driven 
and nozzle gases were modelled using the species N 2 , O 2 , NO, N and 0 and a five equation 
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reaction set. The mole fractions of the species of the driver gas were frozen at the values 
given after combustion, as described above. The specific heat of the (frozen) driver gas 
was taken to follow values from the JANAF Tables. 15 A 250 cell coarse grid was used for 
the driver and driven tubes. In addition, a sliding 800 cell grid was used which followed 
the main shock and the interface down the driven tube an then became strationary at the 
downstream end of the driven tube. The sliding grid allowed much higher resolutions to 
be obteained in the region of the shocks and the interface. There is an abrupt area change 
at the main diaphragm. The nozzle is modelled with a 100 cell grid which extends 2 m 
downstream from the end of the driven tube. For the reference case without driven tube 
contraction, the throat diameter is 7.60 cm. For most cases with driven tube contraction, 
the throat diameter is reduced in proportion to the reduced driven tube diameter at the 
downstream end of the tube. The nozzle diameter at the end of the nozzle grid is 76 cm 
and the “no gradient” boundary condition is used there. At beginning of the calculation, 
the primary diaphragm is suddenly removed. Wfien the main shock reaches the secondary 
diaphragm, the pressure rises steeply; when it reaches 5.0/ x 10 d/cm , the secondary 
diaphragm is removed. 

IV. RESULTS 

A. Time Histories at Nozzle Entrance 

Figure 6 shows pressure time histories at the nozzle entrance for cases 1,0, 2,02, 2,1 and 
2,20. First, we note that use of the driven tube contraction technique can produce about 
a doubling of the nozzle reservoir pressure. Second, the different contraction geometries 
have large differences in their ability to maintain the increased reservoir pressure. With the 
contraction in the downtream half of the driven tube, performance is very poor. Somewhat 
improved performance is obtained if the entire driven tube is converging. By far the best 
performance is obtained if the contraction is in the upstream half of the driven tube. In 
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this case, the nearly flat part of the pressure history extends out to ~13 msec after arrival 
of the main shock. For this best case, the increase in reservoir enthalpy due to the driven 
tube contraction is by a factor of ~1.5. From the pressure variations seen in the first 
~2 msec in Fig. 6, it is clear that the shock tunnel operation is not tailored for the cases 
shown there. Tailoring is computationally expensive, typically requiring ~5 additional runs 
with varying initial driven tube pressure per case and hence, has not been done for most 
cases presented herein. Tailoring of an optimized case for a 2 to 1 diameter driven tube 
contraction will be discussed at a later point. A comparison was then made between the 
nozzle entrance pressure histories for the 2,20 case and the 2,40 case with the contraction 
taking place in the upstream quarter of the driven tube. For brevity, the plots are not 
shown here, but in regard to the ability to maintain the high reservoir pressure the 2,40 
case was found to be much inferior to the 2,20 case. 

Further studies are now presented keeping the basic 2,20 location of the contraction 
section, but considering non-linear diameter tapers. We consider cases 2,20, 2,20tl and 
2,20b. These are linear, trumpet and bell shaped tapers, respectively, and the equations 
of the tapers are given in Sec. IIIA, Eqs. (1), (3) and (2), respectively. Figure 7 shows the 
nozzle entrance pressure histories for these three cases. At this point in the study, we are 
searching for the contraction geometry which produces the most nearly constant nozzle 
reservoir conditions. From Fig. 7, we can see that the 2,20b case is roughly equivalent 
to the 2,20 case in this respect, but that the 2,20tl case is substantially better than the 
2,20 case. Since the trumpet case 2,20tl produced the most nearly constant reservoir 
pressure history, this line of attack was pursued further and two additional trumpet cases 
were investigated. For case 2,20t2 (Eq. (4), Sec. IIIA) the basic shape of the trumpet of 
case 2,20tl was maintained, but the trumpet entrance diameter was increased to equal the 
driver tube diameter, thus eliminating the step area change at the main diaphragm. The 
driven tube diameter at the downstream end of the tube is the same for cases 2,20tl and 
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2 , 20 t 2 (and also for all 2,- cases). For case 2,20t3 the taper equation is given below. 


D ( x v 

— =0.5 + 0.5 1 - — 
D 0 V Ax/ 


( 5 ) 


For cases 2,20tl and 2,20t2, there is an abrupt change in wall slope at the downstream end 
of the taper section. For case 2,20t3 the taper section slope, instead, goes smoothly to zero 
at the downstream end of the section. Case 2,20t3 maintains the step area change at the 
main diaphragm typical of all 2,- cases except the 2,20t2 case. Figure 8 shows the nozzle 
entrance pressure histories for cases 2 , 20 t 1 , 2 , 20 t 2 and 2,20t3. The 2 , 20 t 2 case shows about 
a 15% pressure increase over the 2 , 20 tl case (after the initial pressure spiking). However, 
both the 2 , 20 t 2 and the 2 , 20 t 3 show much more variation of nozzle reservoir pressure 
during the critical first ~13 msec than the 2 , 20 tl case. 

One last non-linear taper case was tested, case 2 , 20 h. The taper equation for this case 
is given below. 


This taper equation is a hyperbolic tangent with minimum slope at the beginning and end 
of the taper section and maximum slope at the middle of the taper section. Compared 
to other cases, for case 2 , 20 h there is a relatively steep wall slope at the middle of the 
taper section and fairly early in the CFD calculation it became apparent that choking was 
occuring near the maximum slope point. Hence, this taper was rejected as unsatisfactory 
afid the solution stopped. Hence, no nozzle entrance time histories were obtained for this 
case. We will discuss choking further in Sec. IVB. From the discussion of this section, we 
conclude that the best contraction geometry of the 2 ,- cases studied is the trumpet case 
2 , 20 tl. 

For the optimized case 2 , 20 tl, the tunnel operation was tailored by varying the initial 
driven tube fill pressure. Calculations were run for pressures of 0.65, 0.75, 0.80, 1.20 and 
1.40 times the original driven tube fill pressure of 8.67 x 10 5 d/cm 2 . Figure 9 shows the 


1 + 


tanh (77 — 2 ) 


tank 2 


( 6 ) 


nozzle entrance pressure histories for 0.65, 0.80, 1.00 and 1.20 times the original driven 
tube fill pressure. The best tailoring is seen to occur for an initial driven tube fill pressure 
of 6.94 x 10 5 d/cm 2 or 0.80 times the value used for most run reported herein. For this 
case, the pressure between 6.2 - 6.6 msec is very nearly equal to that after ~/.4 msec. The 
wiggle in between these times may be a computational artifact due to the limited resolution 
of the shock waves and interface. 

Figure 10 shows the nozzle reservoir pressure and enthalpy histories for the 1,0 and 
the 2,20tl cases (both cases are untailored). We note that for the 1,0 reference case, the 
pressure starts to fall off ~18 msec after initial shock arrival. This is in excellent agreement 
with experimental results from the operation of the Ames 16 Inch Tunnel. 10 As mentioned 
in Sec. IIIB, much of this experimental data was taken at conditions closely matching the 
reference case 1,0 conditions except that the experimental data was taken at somewhat 
lower pressure levels. The enthalpy for the 1,0 case is calculated to begin to fall off ^12 
msec after initial shock arrival. For the 2,20tl case, both the pressure and enthalpy are 
calculated to begin to fall off ~13 msec after initial shock arrival. From Ref. 10, the 
actual driver gas free test time for the cases nearly matching case 1,0 is estimated to be 
~5 msec. This relatively short test time is due to a number of effects which cause the 
driver-driven gas interface to spread out significantly from the ideal sharp interface (see 
Sec. I). The increased pressure and enthalpy of the case 2,20tl nozzle reservoir condition 
are thus maintained sufficiently long to allow full use of this condition up to the likely time 
when driver gas will begin the flow through the nozzle, ending the useful test time. Going 
from the reference case 1,0 to the 2,20t 1 case provides an increase in reservoir pressure 
from 8.01 x 10 8 d/cm 2 to 1.67 x 10 9 d/cm 2 , a factor of 2.09. The reservoir enthalpy is 
increased from 9,500 J/gm to 14,300 J/gm, a factor of 1.51. 

We now discuss cases with a driven tube diameter contraction ratio of 3 to 1. Only 3 
cases were computed. Two of these cases were with linear diameter tapers, cases 3,02 and 
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3,20. The third case was a trumpet case denoted by 3,20tl. The taper of the 3,30tl is of 
the same form as that of the 2,20tl, with constants adjusted to provide a 3 to 1 diameter 
ratio. The taper equation of case 3,30tl is as follows. 


— = 1 - 0.667 
D 0 



( 7 ) 


The last case was chosen because of the superior performance of case 2,20tl. Figures 
11 and 12 show the pressure and enthalpy histories at the nozzle entrance for cases 1,0 
(reference case), 3,02, 3,20 and 3,20tl. The performance of the 3,02 case is very poor, 
considering the 3 to 1 contraction ratio, providing only conditions about equal to those of 
the 2,20 case. The 3,20 case provides much larger gains in enthalpy and pressure. The 
trumpet case 3,20tl provides a further gain in enthalpy of ~10% and a considerably flatter 
pressure history. Choking effects are quite strong for the 3,20 case, but much less serious 
for the 3,20tl case, due to the curved trumpet taper of the latter. The reduction of choking 
effects for the 3,20tl case is likely responsible for its superior performance. (Choking effects 
will be discussed further in Sec. IVB.) It is likely that significantly better taper profiles 
could be found for 3 to 1 diameter contractions; however, limitations of time and available 
computational resources have prevented us from pursuing such cases further at the present 
time. 

Let us examine the performance of the 3,20tl case. From Fig. 11, if we imagine this 
case to be tailored, it seems likely that the nozzle reservoir pressure could be maintained 
constant within ~5% for about 5 msec. Further optimization might well allow nearly 
constant pressure to be maintained for a longer time. Since the arrival of driver gas at the 
nozzle entrance has been estimated to occur ~5 msec after initial shock arrival, it follows 
that full utilization of the tunnel test time capability can be achieved using 3 to 1 diameter 
contraction as well as using 2 to 1 diameter contraction. Going from the reference case 
1,0 to the 3,20tl case provides an increase in reservoir pressure from 8.01 x 10 8 d/cm 2 to 
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2.59 x 10 9 d/cm 2 , a factor of 3.23. The reservoir enthalpy is increased from 9,500 J/gm to 
19,800 J/gm, a factor of 2.08. 

B. Plots of data along driver and driven tubes 

We now turn to plots of data along the driver and driven tubes just before the main 
shock reaches the nozzle entrance. We consider first the the reference case without driven 
tube contraction (1,0), the optimized case for 2 to 1 diameter contraction (2,20tl) and the 
best case for 3 to 1 diameter contraction (3,20tl). This last case uses the trumpet taper of 
case 2,20tl modified for 3 to 1 diameter contraction. These cases all have an initial driven 
tube fill pressure of 8.67 x 10 5 d/cm 2 and are untailored. Figures 13 and 14 shows plots of 
velocity and pressure for these three cases. In these figures, the driver tube extends fiom 
0 to 21 m and the driven tube from 21 to 47 m. From Fig. 13 we see that the gas velocity 
behind the main shock rises from 3.1 to 3.8 to 4.6 km/sec for diameter contraction ratios 
of 1, 2 and 3, respectively. The corresponding pressures behind the main shock rise from 
0.9 x 10 8 to 1.6 x 10 8 to 1.9 x 10 8 d/cm 2 . 

From Fig. 14, even in the downstream one-third of the driver, the pressures are 
somewhat higher (10 - 15%) for the cases with contraction. The pressure drop at the main 
diaphragm station is substantially lower for the cases with contraction. This is due to the 
lower Mach numbers at this station for cases with contraction. The Mach number (plots 
not shown) on the driven tube side of the diaphragm station is ~0.5 for the cases with 
contraction compared to 1.0 for the constant driven tube area case. Further, on a fractional 
basis, the pressure drop in the converging part of the driven tube is much smaller for the 
cases with contraction than for the corresponding part of the constant area driven tube. 
In fact, for the 3 to 1 contraction, there is actually a net pressure rise in the contraction 
section. Finally, for cases with contraction, in the constant area part of the driven tube, 
the pressure falls off, but the performance of the geometries up to this point has been 
sufficiently good that the net pressure behind the main shock still considerably exceeds 
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that for the constant area driven tube case. 

We now discuss choking in the driven tube. Figures 15 and 16 show plots of pressure 
an Mach number along the driver and driven tubes for cases 1,0, 3,20 and 3,20tl. Case 3,20 
is a 3 to 1 diameter contraction case with linear taper and case 3 ,20t 1 is the corresponding 
case with a trumpet taper. For the 1,0 (reference case), the Mach number profile along 
the driven tube shows a smooth increase from 1.0 at the main diaphragm station to ~1.9 
at the end of the expansion fan 3 in Fig. 2 (see Sec. I) and then remains constant out to 
the interface, across which it jumps to ~2.3. In constrast, for case 3,20, severe choking is 
apparent at the end of the driven tube contraction section. There is a sharp pressure peak 
reaching 1.2 x 10 9 d/cm 2 , and Mach number, instead of rising in the contraction section, 
falls from a maximum of ~0.9 to ~0.8 at the end of the contraction section. The rate of 
contraction of the linear taper in the logarithmic sense, d(ln(Area))/ d(x/L), appears to be 
too great at the end of the contraction for a 3 to 1 diameter ratio. (L is the length of the 
driven tube.) By switching to the trumpet taper, case 3,20tl, the choking problem in the 
driven tube is very much ameliorated. The pressure peak at the end of the contraction has 
been reduced by a factor of 2 and the Mach number in the contraction section now rises 
smoothly from ~0.5 to ~1.3. There still does remain a much reduced pressure peak at the 
end of the driven tube contraction section; it is possible that this peak could be removed 
by further optimization of the contraction taper shape. 

.- In Sec. IVA, we noted that case 2,20h, with a hyperbolic tangent taper profile, pro- 
duced severe choking. Figure 17 shows plots of pressure and Mach number along the driver 
and driven tubes for this case. The choking is shown by the sharp pressure peak and Mach 
number drop (from 1.1 to 0.6) near the center of the taper section, where the slope of the 
tube wall is the steepest. For this case also, it appears that d(ln(Area))/d(x/L) is too 
great to avoid substantial choking effects. 

C. Test times and impact pressures in the test section 


16 



First, we consider cases 1,0 (reference case) and 2,20tl (optimized case with 2 to 1 
diameter contraction). As we have noted previously, the nozzle reservoir pressure and 
enthalpy of the 2,20tl case are 2.09 and 1.51 times higher, respectively, than the values 
for the 1,0 case. However, these gains were not achieved without some cost. The volume 
of compressed test gas in the driven tube reservoir is considerably smaller for the 2,20tl 
case than for the 1,0 case. At this time, we have no way to know whether the spreading 
out of the driven-driver gas interface (see Sec. I) will be better, worse or about the same 
for the 2,20tl case, compared to the 1,0 case. Lacking this information, we shall assume 
that the volume of compressed driver-gas-free test gas at the nozzle entrance for the 2,20tl 
case will be the same fraction of the driven tube volume as for the 1,0 case. The volume 
of the driven tube for the 2,20tl case can readily be calculated, using the taper equation 
(3) to be 0.354 of the volume for the 1,0 case. 

We now introduce a number of ratios, all of which are formed by taking the parameter 
for the case with driven tube diameter contraction (here, case 2,20tl) and dividing by the 
same parameter for the reference case (1,0). R p is the ratio of driven tube reservoir 
pressures, here equal to 2.09. R/j is the ratio of driven tube reservoir enthapies, here equal 
to 1.51. R„ is the ratio of driven tube reservoir gas volumes, here equal to 0.354. R t , R a th, 
Rata, R dts and R p i are the ratios of test time, nozzle throat area, test section area, test 
section diameter and test section impact pressure, respectively. Now, R p , R/ t and R„ are 
given for the case 2,20tl - case 1,0 comparison. However, Rath and Rais can be varied and 
allow one to trade off among R a ts, Ri and R pi- R<> R»> Rath and R h can be related as 


follows. 


7 ? = — 1 

Rath R'h 


Equation (8) can be obtained by noting that the flow time of the test gas slug through 
the nozzle throat is proportion to the volume of the gas slug and inversely proportional 
to the area of the nozzle throat and the sound speed of the gas. Further, the following 
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approximate relation for hypersonic flow can be written. 


R 


V' 


Rath t-) 

~R ™ p 
-K-ats 


( 9 ) 


Equation (9) can be obtained, for example, from the discussion of the treatment of the 
impulse function in Ref. 16 by assuming that the Mach number is much greater than unity. 
Combining Eqs. (9) and (10), we obtain 


RtRjnRats = = °' 603 


( 10 ) 


The numerical constant in Eq. (10) was obtained using R^, Rp and R/i for the case 2 ,20 1 1 
- case 1,0 comparison (in addition to one case for 3 to 1 diameter contraction). 

Using Eqs. (9) and (10), a number of cases will be constructed and discussed. The 
parameters for these cases are shown in Table 2. 


Table 2 Test times and impact pressures 


Nozzle 

Driven 

R ath 

Ra^ 

R dts 

R/i 

Rt 

R pt 

case 

tube 







# 

case 







i 

1,0 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

2 

2,20tl 

0.25 

0.25 

0.50 

1.51 

1.15 

2.09 

3 

2,20tl 

1.00 

1.00 

1.00 

1.51 

0.289 

2.09 

4 

2,20tl 

0.50 

1.00 

1.00 

1.51 

0.577 

1.04 

5 

2,20tl 

0.25 

1.00 

1.00 

1.51 

1.15 

0.522 

6 

2,20t 1 

0.50 

0.50 

0.707 

1.51 

0.577 

2.09 

■- 7 

3,20tl 

0.367 

0.367 

0.606 

2.08 

0.420 

3.23 


Nozzle case 1 is simply the reference case without driven tube contraction. All ratios for 
this case are, by definition, equal to unity. For the remaining nozzle cases, driven tube 
contraction, cases 2,20tl or 3,20tl is assumed. For nozzle case 2, the nozzle throat and test 
section linear dimensions are half the values for the reference case. The test section impact 
pressure is increased by a factor of ^2 and the test time is actually increased 15% above 
the value for case 1,0. The disadvantage is that the model must be twice as small as for the 
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reference ca.se. For nozzle case 3, the nozzle throat and test section linear dimensions kept 
equal to the values for the reference case. The test section impact pressure is increased by 
a factor of two and the model size is kept undiminished, but the test time is ~0.29 that for 
the reference case. For many cases, this would be unsatisfactory, since there might well not 
be sufficient time for the flow the stabilize properly in the nozzle and/or over the model. 
Nozzle cases 4 and 5 are cases where the nozzle throat area has been reduced by factors of 
2 and 4, respectively, from that of case 3, while keeping the test section size undiminished. 
The result is that the test times are increased in inverse proportion to the throat areas, 
but the impact pressures in the test section are correspondingly reduced. 

For nozzle cases 2 - 5, we have not considered that less test time should be required 
to establish the flow if the nozzle and model are smaller and also that less test should be 
required when the free stream velocity is faster due to increased reservoir enthalpy. To 
allow for this we now assume the following. 

R'L = RhR: ( u ) 

The condition of Eq. (11) is that the test section (and presumably nozzle and model) 
dimensions should vary, for the various cases studied, proportionally to the free stieam 
velocity and the test time. Applying the condition of Eq. (11) to Eq. (10) yields the 
following equation. 

R pi R 3 t = 0.401 (12) 

Let us assume that we wish to obtain the full pressure advantage of the driven tube 
contraction technique and therefore take R pt = 2.09. From Eqs. (9), (11) and (12), we 
then obtain R ( = 0.577, R a( * = R a th = 0.500 and D a<s = 0.707, nozzle case 6 in Table 2. 
For this case, the test time is 0.577 times that of the baseline case, but this shorter time 
should be just as effective a setting up the flow as in the reference case, since the model 
is smaller and the flow is faster. The test section and model size are 0.707 that for the 
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reference case and the full pressure advantage of the driven tube contraction technique is 
available in the test section. 

We should note here that nozzle cases 4 and 5 in Table 2 involve a change in area 
ratio of the nozzle, since R at3 is not equal to R at h for these cases. If these cases were to be 
used, there would be a tendency for the test section Mach number to change from that for 
the reference case. For example, assuming a reference case test section Mach number of 
5 and and effective 7 value of 1.25, doubling and quadrupling the nozzle area ratio would 
increase the Mach number to ~5.6 and ~6.2, respectively. 

We now discuss the case for a 3 to 1 diameter contraction, using cases 3 , 20 tl and 
1,0. R p , Rh and R v for this case are 3.23, 2.08 and 0.222, respectively. For brevity, we 
go directly to the case where we maintain the full pressure advantage of the contraction 
technique and allow the test time to be shortened corresponding to a shorter nozzle and 
model and faster flow. The resulting parameters are shown in Table 2, nozzle case 7. The 
test time is 0.420 times that of the baseline case but should be equally satisfactory at 
setting up the flow properly. The model dimensions are 0.606 times those for the baseline 
case. Other possible combinations of R a js and R a th could of course be used, allowing a 
wide variety of cases to be obtained, including those corresponding to nozzle cases 2 - 5 
for 2 to 1 diameter contraction. 

D. Further work 

Pursuing the converging driven tube technique further would likely involve the follow- 
ing steps. First, a contraction ratio yielding approximately the desired condition would be 
selected. Then, inviscid analyses similar to those presented here should be done, but in 
more depth, to allow the best location and shape of the contraction section to be found. 
Then, axisymmetric viscous flow analyses with turbulence modelling should be carried out. 
This would allow possible additional choking effects due to boundary layers to be studied 
as well as some of the effects which cause spreading out of the driver-driven interface. 
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The latter include jetting of the driver gas upon initial rupture of the main diaphragm, 
effects of the boundary layer upon the motion of the initial shock and the driver-driven gas 
interface, bifurcation of the reflected shock in the boundary layer and, possibly, instability 
effects at the interface. The best viscous flow cases should then be tested experimentally 
in a small facility. Finally, if the predicted gains in pressure and enthalpy are verified 
experimentally at small scale, one might consider the modification of a major facility to 
use the driven tube contraction technique. 

V. CONCLUSIONS 

High stagnation pressures and enthalpies are required for the testing of aerospace ve- 
hicles such as aerospace planes, aeroassist vehicles and reentry vehicles. Shock tunnels are 
among the most useful facilities to perform such tests. For large, high pressure production 
facilities, the most practical ways of achieving high enthalpy in the driver of a shock tunnel 
are to use light gases, hydrogen or helium, and to heat the gas using electrical resistance 
heaters, combustion or piston compression. The enthalpies achieved in these ways are lim- 
ited and there are also pressure limits on the driver tube. Given a fixed driver condition, 
the nozzle reservoir enthalpy and pressure in the driven can be varied by changing the 
size of the driven tube and the initial driven tube fill pressure. The gains in pressure and 
enthalpy achievable using the former technique are very limited. The latter technique can 
be used to achieve significant increases in enthalpy, but at the cost of decreased pressure 
and test time. 

A new technique has been presented which allows substantial increases in nozzle reser- 
voir pressure and enthalpy to be achieved simulataneously. In this technique, part of 
the driven tube is made converging. This technique has been investigated using a one- 
dimensional, inviscid CFD code with chemical kinetics. The computational techniques 
were applied to analyze possible modifications of the NASA Ames 16 Inch Hypersonic 
Shock Tunnel. Contraction ratios for the driven tube diameter of 2 to 1 and 3 to 1 were 
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studied, the former in more detail. The length, location and shape of the converging sec- 
tion have been varied and partially optimized cases, within the limitations of available 
computer time, have been found. Data has been presented in the form of pressure and 
enthalpy histories at the driven tube nozzle reservoir and as snapshots of velocity, pressure 
and Mach number along the the driver and driven tubes. 

For a driven tube diameter contraction ratio of 2 to 1, increases in nozzle reservoir 
pressure and enthalpy of 2.1 and 1.5, respectively, were obtained for the best cases studied. 
For a driven tube diameter contraction ratio of 3 to 1, the corresponding numbers were 
found to be 3.2 and 2.1. Choking in the driven tube was found to be a problem for some 
cases studied. It was found that choking effects could be minimized by the proper selection 
of the location, length and shape of the contraction section. Using both 2 to 1 and 3 to 1 
diameter contraction ratios, test times during which the nozzle reservoir maintained nearlv 
constant pressure and enthalpy were evaluated. These were found to exceed or equal the 
experimentally estimated driver-gas-free test time for the Ames Shock Tunnel. Hence, the 
f ull test time of the tunnel would appear to remain usable if it were modifed with optimized 
2 to 1 or 3 to 1 driven tube contraction profiles. 

The gain in nozzle reservoir pressure and enthalpy achieved using the driven tube 
contraction technique is not obtained without cost. There is a significant reduction in the 
available volume of gas at the nozzle reservoir conditions. For this reason, there will a 
tendency for the test time and the size of the test section which can be operated at a given 
impact pressure to be reduced when the driven tube contraction technique is used. By 
varying the nozzle throat and test section dimensions, trade-offs can be made between test 
time, test section size and test section impact pressure. These questions have been studied 
for the best 2 to 1 and 3 to 1 contraction cases investigated. For these cases, for one set of 
calculations, it was assumed that when using the contraction technique, it was desired to 
maintain the full gain in impact pressure in the tunnel test section. Also, allowances were 
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made for the fact that the test section flows could be stabilized in shorter times since the 
models and nozzle were smaller and the flows faster (on account of the higher enthalpies) 
when contraction was used. For these conditions, for the best 2 to 1 contraction case, 
the model size was found to be 0.707 of that for the baseline case (without contraction) 
and the test time was found to be 0.577 of that for the baseline case. For the best 3 to 1 
contraction case, the corresponding numbers were found to be 0.606 (for model size) and 
0.42 (for test time). However, we must point out again that these shorter test times should 
be fully sufflcent to stabilize the flows as well as in the baseline case without contraction. 
Other combinations of nozzle throat area and test section area have also been discussed. 
For example, for the best 2 to 1 contraction case, the full impact pressure gain can be 
maintained and the full model size maintained, with the penalty that the test time drops 
to 0.29 of that without contraction. 

APPENDIX A 

This appendix describes a simple zoned, inviscid, ideal gas method of estimating the 
performance of the driver and driven tubes of a shock tunnel. We refer to the X-T diagram 
of Fig. 2 of the main text. A description of this diagram is given in Sec. I of the main text 
and will not be repeated here. The gases are taken to be ideal, with values estimated from 
the enthalpy values of the JANAF tables for roughly estimated temperature ranges for the 
various zones. With one exception, all friction and heat transfer effects are neglected. The 
driver gas is initially taken to be 2H 2 + 0 2 + 9He, which burns to 2H 2 0 + 9He. The 
driven gas is air, taken as 3.77N 2 + 0 2 . The one dissipative effect which is modelled is 
that 30% of the energy which would be released upon complete driver gas combustion is 
assumed to be lost. This was found to bring the calculated shock velocities into rough 
agreement with experimental values and thus crudely allows for dissociation, heat transfer 
and frictional losses. Expansion waves, such as El and E2 and W (if it is an expansion 
wave), are treated as simple centered expansion wave systems. The flow from 3’ to 3” 
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is taken to be a simple steady expansion to Mach 1 at 3 . Shocks Si, S2, the shock 
above and to the left of M and W (if it is a shock) are treated by standard shock wave 
analyses. The flow velocity in zones 2 and 3 is varied until the pressures in these two 
zones match, to obtain the correct velocity. A similar procedure is applied in the two 
zones above M to obtain the correct velocity and pressure of these zones. To tailor the 
operation of the tunnel, the initial pressure in zone 1 is varied until wave W is of zero 
strength. For equilibrium interface operation, the conditions in zone 5 are first calculated. 
Then the pressure which zone 3 would reach if brought to rest by a single shock (p3, shoe/.-) 
is calculated. The equilibrium interface nozzle reservoir conditions are then calculated 
by assuming isentropic compression or expansion from state 5 to the pressure ps^hock- 
The code available at Ames can allow for flow through a shock tunnel nozzle of specified 
area. The calculation is somewhat shorter if the nozzle is taken to be plugged. Results for 
plugged nozzles were used in the discussion of the main text. 
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Figure Captions 


Fig. 1 Sketch of shock tunnel (not to scale). 

Fig. 2 x-t wave diagr am for shock tunnel. SI and S2 denote shock waves, El and 
E2 denote expansion wave systems, HE denotes the head of the expansion wave system 
reflected from the end of the driver, I denotes driver-driven gas interface. 1,2, 3, 3’, 3” ,4, etc. 
denote thermodynamic states. Also, see text. 

' Fig. 3 Effect of driven/driver area ratio (A<*„/A dr ) on ratio of reservoir pressure to 
driver pressure (Prea/Pdr)- Tailored interface. See text for other details of operating 
conditions. 

Fig. 4 Non- tailored operation of shock tunnel. Abscissa is ratio of driven tube fill 
pressure to that required for tailored interface operation. On the ordinate are shown the 
ratios of nozzle reservoir pressure (R p ) and enthalpy (R*) to those obtained for tailored 
interface operation. R r is the volumetric compression ratio of the doubly shocked gas in 
the nozzle reservoir divided by the corresponding value for tailored interface operation. 
See text for details of operating conditions. 

, Fig. 5 Various driven tube contraction profils studied herein. Diameters are to con- 
sistent scale, length is not at same scale as diameters. 

Fig. 6 Pressure histories at nozzle reservoir for cases 1,0, 2,02, 2,1 and 2,20. 

Fig. 7 Pressure histories at nozzle reservoir for cases 2,20, 2,20b, and 2,20tl. 

Fig. 8 Pressure histories at nozzle reservoir for cases 2,20tl, 2,20t2 and 2,20t3. 


Fig. 9 Pressure histories at nozzle reservoir for case 2,20tl. Parameters shown on 
curves axe ratios of driven tube initial fill pressure to that for the first case calculated. The 
tailoring effect is shown. 

Fig. 10 Pressure and enthalpy histories at nozzle reservoir for cases 1,0, and 2,02tl. 

Fig. 11 Pressure histories at nozzle reservoir for cases 1,0, 3,02, 3,20 and 3,20tl. 

Fig. 12 Enthalpy histories at nozzle reservoir for cases 1,0, 3,02, 3,20 and 3,20tl. 

- Fig. 13 Snapshot along driver and driven at time initial shock reaches nozzle entrance. 

Velocity, cases 1,0, 2,20tl and 3,20tl. 

Fig. 14 Snapshot along driver and driven at time initial shock reaches nozzle entrance. 
Pressure, cases 1,0, 2,20tl and 3,20tl. 

Fig. 15 Snapshot along driver and driven at time initial shock reaches nozzle entrance. 
Pressure, cases 1,0, 3,20 and 3,20tl. 

Fig. 16 Snapshot along driver and driven at time initial shock reaches nozzle entrance. 
Mach number, cases 1,0, 3,20 and 3,20tl. 

Fig. 17 Snapshot along driver and driven at time initial shock reaches nozzle entrance. 
Pressure and Mach number, case 2,20h. 
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ABSTRACT 

A series of hydrogen-air mixing and combustion 
experiments were conducted at high enthalpy 
conditions at the NASA-Ames Research Center 
16-Inch Shock Tunnel. The goals of the tests 
were to demonstrate the facility capability for 
high speed propulsion testing and to obtain 
limited data on the performance of 30° flush 
wall injectors. The experimental results were 
compared with CFD analysis. 

INTRODUCTION 

The rebirth of hypersonics research in the U.S. 
has stimulated the need for ground test 
simulation of supersonic combustion ramjets at 
high flight Mach numbers. Until recently, very 
little data was available at Mach numbers 
greater than about 8. Ground testing above 
Mach 8 requires high pressure, high enthalpy 
facilities such as shock tunnels^' 3 , and 
expansion tubes. ^ These facilities have test 
times which are free of driver gas 
contamination on the order of a millisecond. A 
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notable exception to the short duration, high 
enthalpy facility tests is a recent set of 
experiments in the Mach 9-12 regime at the 1 00 
MW arcjet facility at NASA-Ames 3 . 

As the need for data at high enthalpy 
conditions becomes more important in the design 
of scramjet engines, additional facilities are 
being brought on line for propulsion testing. 
This paper presents the results of the first 
high speed propulsion related experiments in 
the Ames 16-Inch Shock Tunnel simulating 
combustor inlet conditions at approximately 
Mach 14. The test program was designed as an 
initial entry to demonstrate the facility 
capability and to obtain limited data on the 
characteristics of 30° flush wall hydrogen fuel 
injectors. 

FACILITY 

A brief description of the Ames 16-Inch Shock 
will be presented. Complete details of its 
operation and flow characteristics can be found 
in Refs. 6 & 7. The NASA-Ames 16-Inch Shock 
Tunnel (Fig. 1) has a 21 m long, 43 cm diameter 
driver tube and a 26 m long 30 cm driven tube. 
The facility receives its name from the 16 inch 
naval cannons that were used to construct the 
driver tube. The nozzle is contoured and 
designed for a nominal Mach number of 7 at low 
enthalpy, however due to the high enthalpy 
conditions the actual free stream Mach number 
was about 5.8. The facility was originally built 
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as counterflow to a ballistic range and not as a 
stand alone shock tunnel. The test section in 
this entry was at the last station of the 
ballistic range and was located 2.6 m 
downstream of the nozzle exit. The test section 
diameter was 0.99 m. For the present tests, the 
nozzle was operated with a test section to 
throat area ratio of 170; however, the nozzle 
can be run at area ratios as low as 95. 

The combustion driver is operated with a 2H2 + 
1.1702 + 9.36He gas mixture. The gas is ignited 
by using four 0.019 cm diameter tungsten wires 
strung the length of the driver tube. The wires 
are connected in parallel and heated to 2000 K 
by the discharge from a capacitor bank. The 
ignition produces smooth combustion with a rise 
time of approximately 25 ms and produces a 
pressure rise of approximately 8:1. The 
maximum pressure ratings of the driver and 
driven tubes are 680 atm. For the data 
presented here, the final after bum pressure in 
the driver tube was 272 atm. Flat, scored 
diaphragms are used which self break at the 
completion of the driver bum. The test time 
which is free of driver gas contamination for 
the conditions used here was estimated to be 3-5 
ms. 
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The tests were performed on a wedge model 38 
cm long, with a 46 cm wide leading edge and 
with the instrumented surface inclined 11° to 
the free stream flow. The model is shown in 
Fig. 2. A spanwise row of five sonic 0.305 cm 
diameter injector ports was located at 22 cm 
from the leading edge. The injectors were 
inclined at 30° to the model surface and were 
8.33 diameters apart. 


The gages were arranged mainly in three 
streamwise rows with a 3.81 cm streamwise 
separation. The center row of transducers and 
the trailing edge pitot probe were in line with 
an injector port and the outboard rows were 
midway between ports. 

Data Acquisition 

The data-acquisition system (DAS) included 68 
channels of 1-MHz digitizers (LeCroy 8210) 
with 16 kbytes memory per channel, and 20 
channels of digitizers (LeCroy 6810) with up to 
5-MHz sampling rate and 512 kbytes of memory 
per channel. The system was controlled by a 
stand-alone Macintosh II fx computer. 
Amplifiers with gains of 10 and 100 and 
frequency responses of 70 kHz were used. 

Flow Visualization 

A versatile laser holographic interferometer 
(see Fig. 3) was used for flow visualization in 
the spanwise direction. The technique used was 
a dual-plate, double-pass system. The dual- 
plate technique provides system insensitivity 
to optical nonuniformities and provides 
flexibility in reconstruction. The double-pass 
technique provides increased sensitivity at 
lower pressures. The light source was an 
injection-seeded, frequency-doubled, Nd:YAG 
laser producing pulses at 10 Hz, 532-nm 
wavelength, 200-mJ pulse energy, and 6-nsec 
pulse width. The beam distribution optics were 
located directly adjacent to the optical access 
port at the test section, while the laser light 
source was remotely located and directed to the 
test section by a series of mirrors. 

TRANSVERSE TNTECTION 


The model instrumentation consisted of pitot 
probes, static pressure gages and heat flux 
gages. At the leading edge there was a two 
probe pitot rake and at the trailing edge there 
was a seven probe rake. A total of 12 Kulite 
surface pressure transducers were used, as shown 
in Fig. 2. In addition there were two hydrogen 
injector plenum pressure transducers, and two 
hydrogen flow venturi pressure transducers. 
The heat transfer gages included 10 Type E 
coaxial thermocouples, 2 thin skin slug 
calorimeters and two platinum thin film gages. 


Transverse injection of hydrogen through 
discrete holes for application in scramjet 
combustors has been investigated by several 
authors, mostly in low enthalpy continuous 
flow facilities 8 ' 12 . Due to the long run time in 
these types of facilities, in-stream 
measurements such as gas sampling, pitot 
surveys and static pressure surveys were 
possible. These results have been used for 
extensive calibration of CFD codes. Since data 
at high enthalpy conditions is very sparse and 
in-stream data such as species concentrations is 
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very difficult to obtain, it would be desirable if 
the cold flow mixing results could be extended 
to high enthalpy conditions. In one of the low 
enthalpy transverse injection experiments, 
McClinton 10 investigated an array of five 
injectors inclined at angles from 30° to 90° and 
found the best mixing and penetration for 
lowest injection angles. The lower angles were 
also found to produce less total pressure loss. 
Another significant advantage of injection at 
low angles for scramjet applications is that the 
downstream axial momentum of the jets can 
produce a sizable fraction of the net engine 
thrust. Thus, injection at 30° was chosen for 
this investigation. 

TEST CQM9mQN5 jkMAXRIX 

Tests were performed at the test conditions 
shown in Table 1. The shock tube conditions 
were calculated using standard equilibrium 
shock tube relations. The calculation used the 
measured shock speed and driven tube initial 
conditions to compute the conditions behind the 
reflected shock, and then performed an 
isentropic expansion from the calculated 
pressure to the measured pressure to account for 
the non-tailored interface. The enthalpy of 9.5 
MJ/kg is roughly equivalent to the total 
enthalpy at a flight Mach number of 14. The 
nozzle exit conditions were calculated using a 
quasi-l-dimensional nonequilibrium code. 13 
The nozzle exit Mach number was 
approximately 5.8 for these conditions and the 
flow behind the shock wave had a Mach 
number of 4.4. The static pressure on the surface 
was 0.25 atm and the temperature after 
compression was about 1800 K. The boundary 
layer thickness was estimated to be on the 
order of 1 mm. Though the pressure is lower 
than typical, the velocity, temperature and 
Mach number are close to the combustor inlet 
conditions of a scramjet engine at Mach 14. 

Three types of tests were run: runs with air flow 
and no injection (tare); runs with H 2 injection 
into N 2 flow (mixing), and runs with H 2 
injection into air (combustion). An important 
parameter regarding the penetration of 
transverse jets is the momentum flux ratio of the 

injectant to that of the free stream, i.e., q 
=(pV 2 I/pV 2 0 o) / where V is the velocity, and p 
is the density. The fuel will penetrate further 


into the free stream for higher values of this 
ratio. For these experiments, the dynamic 
pressure of the hydrogen jets was 
approximately equal to that of the shocked 

flow parallel to the wedge (q = 1 condition) for 
all injection runs except one, for which it was 

twice as high (q = 2 condition). The fuel total 

pressure for the nominal q = 1 condition was 8.8 
atm and the discharge coefficient of the 
orifices was 0.78. Ten runs, including several 
repeated runs, were made. 


Code Description 

The NASA Langley Research Center (LaRC) 
SPARK series of codes were used to perform the 
computational analysis. The two-dimensional 
elliptic, explicit, finite difference code was 
originally developed by Drummond (Ref. 15). 
This code was subsequently extended to three 
dimensions by Carpenter (Ref. 16) and this has 
been converted into a parabolized version by 
Kamath (Ref. 17). The codes have been 
extensively validated for a variety of flow 
types. Specifically, the codes have been used 
to model experiments involving hydrogen 
injection into a supersonic stream (Refs. 18-20). 
The codes have a choice of a second order 
MacCormack, a fourth order cross MacCormack, 
or a fourth order Gottlieb MacCormack 
algorithm to solve the mass, momentum, 
energy, and individual species mass 
conservation equations. The cross MacCormack 
algorithm was used exclusively in this 
investigation. This algorithm achieves fourth- 
order spatial accuracy at steady state 
although it is only second-order accurate in 
space and time during the transient 
development. Additionally, fourth order 
damping terms were used which are based on 
pressure and temperature gradients. The 
elliptic codes can be run in a time accurate or a 
local time stepping mode which accelerates the 
convergence to a steady state solution. 
Unfortunately, at this time it is impractical to 
solve a three-dimensional problem with a 
"real time" approach. Also, the programs 
include an internal grid generation capability 
developed by Smith & Weigel (Ref. 21) which 
were used in this analysis. 
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The SPARK codes employ a variety of 
turbulence models. For the analysis contained 
in this paper, the Baldwin-Lomax (Ref. 22) 
algebraic model was used. This method has 
been widely used due to its ease in application. 
Algebraic methods simply model turbulence as 
an increase in the transport coefficients. To 
determine the turbulent portion of the 
viscosity, or eddy viscosity, a length and 
velocity scale must be determined. This model 
bases the velocity scale on the vorticity 
distribution and the length scale on the 
distance from the wall. Although this model is 
not appropriate for injection flowfields, it has 
been calibrated to some degree for jet mixing 
through the use of the turbulent Schmidt 
number as reported by Riggins (Ref. 23). 
Recently, studies have been performed 
successfully using a turbulent Schmidt number 
of 05 (Refs. 18-20) and this value has been 
utilized in this study. Also, to prevent the jet 
induced vorticity from creating unphysically 
high turbulent viscosities the eddy viscosity 
was limited to 1000 times the laminar 
viscosity. Other models such as the k-e and 
Reynolds stress have recently been developed 
and are now in place in some versions of the 
code. Studies are now underway to evaluate 
these models for injection and combustion 
flowfields. 

The SPARK codes have a generalized 
chemistry package wherein the source term can 
be treated implicitly and the capability exists 
to include any number of reactions The reacting 
portion of this work was performed using the 
seven-reaction-seven-species finite rate model 
which is a subset of a model used by Drummond 
(Ref. 15). The reaction rate constants for the 
chemistry model are given in Table 2. For all 
the- calculations, SPARK carries the nitrogen as 
an inert species and its value is determined by 
subtracting the summation of the mass fractions 
of all other species from one.. 

Computational Approach 

The computational approach taken in this 
work was typical of injection flowfields. It 
consists of a sequence of axial blocks as 
illustrated in Fig. 4. The purpose of the 
blocking is to only use a level of code 
complexity that is required in any given region 
of the flow. The inflow conditions are as 


specified in the figure. The free stream 
velocity was broken into two components, one 
parallel and one perpendicular to the wedge 
surface to facilitate the modeling of the wedge. 
For the first zone, the PNS code was run from 
the leading edge of the test section to five 
injector diameters upstream of the jet. Since the 
geometry has no lateral variation upstream of 
the injection point, this region was solved two- 
dimensionally. The second zone was solved 
using the three-dimensional elliptic code. The 
injection conditions are listed on Fig. 4. All of 
the elliptic cases were solved using local time 
stepping to accelerate the convergence to a 
steady state. As in many problems of this type, 
the character of the reacting flowfield is no 
better than quasi-steady so convergence 
requirements for these cases are in reality 
relatively unchanging mean values for wall 
functions such as pressure and skin friction and 
cross-sectional integrated parameters such as 
fuel mixing efficiency and mass conservation. 
The outflow plane from the elliptic solution 
was then passed to the final zone at a plane 
five injector diameters downstream of the 
injection point. Finally, the 3-D PNS code was 
run from this point to the end of the wedge 

Grids and Boundary Conditions 

For the upstream PNS solution, a grid of 5 x 80 
was used. The domain height was set large 
enough to completely capture the wedge 
leading edge shock. These grids were generated 
using the capability in SPARK and were 
clustered along the wedge surface. The wedge 
surface was treated as a no-slip boundary, with 
a wall temperature of 300 K and non-catalytic 
kinetics. 

The grid for the injection region was of the 
dimension (71 x 51 x 80). The jet orifice was 
modeled as a rectangle using two hundred and 
seventy three nodes with specified 
temperature, pressure, and velocity. The 
discharge coefficient was accounted for by a 
slight change in incoming jet pressure. No slip 
boundary conditions were again applied on the 
wedge surface. This grid uses the outflow from 
the PNS code as an inflow condition. In 
addition, a supersonic boundary condition was 
imposed on the outflow plane The lateral 
domain extended from the jet centerline to 
between centerlines with both lateral 
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boundaries modeled as symmetry planes. The 
resulting lateral domain was 1.27 cm wide. The 
grid used clustering at the wedge surface in a 
vertical sense and at the injector in an axial 
sense. The downstream grid utilizes the same 
cross-section as the elliptic zone and continues 
to the end of the wedge. This grid used the 
plane from the elliptic solution as an inflow 
condition and was clustered in a similar manner 
to the other grids. 

RESULTS 

Facility Performance and Data Quality 

One of the primary goals of the experiments 
was to verify the the facility compatibility for 
high speed propulsion testing. The key issues 
in this regard are flow quality, flow 
contamination, test time and repeatability. 
Typical data traces are shown in Fig. 5 in this 
regard. Fig. 5a is a trace of the nozzle reservoir 
pressure. The trace shows that the pressure is 
roughly constant for a period of 18 ms after 
which the pressure drops due to the arrival of 
the reflection of the expansion fan from the far 
end of the driver tube. However, it is well 
known that mixing of the driver gas and the 
driven gas at the contact surface reduces the 
period of uncontaminated test time. Previous 
results have estimated this period to be from 3- 
5 ms for these conditions in this facility. Fig. 5b 
shows a typical facility nozzle static pressure 
trace located at an area ratio (A/A*) of 45 
showing a region of steady pressure followed by 
a slight drop off. This drop off in static 
pressure tends to correlate with the arrival of 
the driver gas due to the fact that the driver 
gas has a much higher value of the ratio of 
specific heats. Fig. 5c is a plot of the static 
pressure on the model surface which also shows 
a steady period of around 5 ms. The pitot 
pressure trace in Fig. 5d, taken at the model 
trailing edge, also shows a reasonably constant 
period of 5 ms. At most transducer locations, 
variations of static and pitot pressures and 
heat flux between the repeat runs ranged from 
2% to 7%. 

Reference 14 summarizes certain rules of thumb 
regarding flow establishment time for pulse 
facility testing. The important parameter is 
the ratio of the slug length, which is the 
expanded extent of the slug of uncontaminated 


gas in the test time t, to the characteristic 
length of the model L. The criteria are tu/L >2 
for attached turbulent flow, and tu/L > 3 for 
attached laminar flow, where u is the flow 
velocity. For the conditions here, the slug 
length based on a test time of 3 ms (which 
represents a conservative estimate of the test 
time ) is about 12 m. Based on this criteria, the 
facility should be able to test large scale 
combustors on the order of several meters long. 
For this experiment, with a model length scale 
of -0.5 m, this criteria is easily met and the 
flow can be said to be fully established in the 
duration of the uncontaminated test time. 

Flow Visualization 

The holographic interferograms were 
reconstructed as shadowgraphs and finite 
fringe interferograms. Examples of finite fringe 
interferograms are shown in Figs. 6 and 7 for the 
cases of tests with and without injection of 
hydrogen into air. Fig. 6 shows that the shock 
from the leading edge of the model falls 
between the fifth and sixth pitot heads, which 
is what was expected from the free stream 
conditions and the turning angle of 11°. The 
fringes are curved in the shock layer most 
likely due to the three dimensional nature of 
the flow field around the model. 

Fig. 7 shows the wedge shock, the boundary 
layer, the hydrogen plume and the shock 
system attached to the plume for the case of 

hydrogen injection into air with q=l. The 
plume shocks begin to merge with the wedge 
shock upstream of the model trailing edge and 
force the wedge shock somewhat upwards for 
cases with hydrogen injection. Fig. 8 is a 
comparison of a shadowgraph reconstruction 
corresponding to the same run as in Fig. 7, with 
a plot of the path averaged density contours 
from the computation. The fuel penetration 
height in the computation is defined as the .005 
H2 mass fraction contour, while the fuel mixing 
layer shown in the shadowgraph corresponds to 
a region of high density gradients. The 
comparison shows that the shock structure in 
the computation agrees with the experiment, 
however it appears as if the fuel penetrates 
further in the computation than the 
experiment. However, this may not be a valid 
comparison since the low concentrations 
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associated with the 0.005 H2 mass contour may 
not be discemable in the shadowgraph. For the 
case of injection with matched dynamic 
pressure, the distance to the edge of the mixing 
layer in the shadowgraph at the pitot rake, h, 
was seen to be about 1.3 cm, which corresponds 
to h/d = 4.3, where d is the jet diameter. The 
penetration of hydrogen plume was seen to be 
substantially greater for the higher injector 

dynamic pressure condition (q=2), as expected. 
The edge of the mixing layer at the pitot rake 
was seen to extend past the third probe, 
approximately 2.2 cm from the surface which 
corresponds to h/ d=7.2. The strength of the 
injector bow shock was also greater than for the 

case with q =1. 

Pitot Pressure 

Figure 9 shows the trailing edge pitot pressure 
profiles for the tare (no injection), mixing and 

combustion runs at q=l and combustion 

conditions at q=2. There were seven pitot 
heads, numbered outwards from the wedge 
surface. The sixth head was inoperative for 
the duration of the tests. All pressure data was 
averaged over 4 ms, starting shortly after 
arrival of the starting shock. The tare run 
shows the wedge shock located outside of the 
fifth probe head, in agreement with the flow 
visualization. The impact pressure remains 
high all the way down to the first probe head. 
For the mixing and combustion runs, the low 
dynamic pressure region of the hydrogen plume 
is apparent adjacent to the wedge surface. It 

covers the lowest two probe heads for q=l and 

the lowest 3 heads for the q=2 condition. The 
plume pitot pressure appears to be about 10% 

lower for the q=l combustion run than for the 
mixing run. This can be explained by the 
higher total pressure loss associated with 
combustion. 

For all runs with injection, the impact pressure 
outside the H2 plume and inside the wedge 
shock is greater than the corresponding values 
from the tare run. This is most likely due to the 
additional compression of the free stream by 
the plume shock system, which is of course not 
present in the tare run. 


Heat Flux Measurements 

Three types of heat flux gages were installed on 
the model for comparison and evaluation 
purposes. Platinum thin film gages are 
commonly used in pulse facilities, however 
previous experience in the facility nozzle 
indicated that these gages might not survive 
the 18 ms period of high impact pressure, 
which is longer than other typical pulse 
facilities. This proved to be true for this 
application also. After the second shot, the 
platinum thin film gages had been destroyed, 
showing infinite resistance. For each run, 
typically less than 2 of the coaxial gages would 
show anomalous behavior. They would give 
good data at the start of the run and then open 
up to infinite resistance. However, these gages 
could be restored to service by abrading them 
with an abrasive eraser and burnishing them 
with a rounded steel rod. At the end of the test 
series, all coaxial gages were still operational. 
One slug calorimeter failed during the ninth 
test run; the other remained operational 
throughout the test series. Based on this 
performance, the coaxial gages were a more 
robust choice for heat flux measurements. 

The two slug calorimeters showed agreement 
with the coaxial thermocouple data to within 
a few percent, however the slug calorimeters 
have a much slower response time (-.5 ms) than 
the coaxial thermocouples, hence only the 
thermocouple data were presented in Figure 10. 
Since coaxial thermocouples gages measure 
surface temperature, the data was reduced to 
heat flux by the standard equations for 
transient heat flux to a semi-infinite slab. An 
automated data reduction scheme was 
incorporated into the data acquisition system 
to reduce the raw data to heat flux. 

Figure 10 shows the heat flux for the tare (no 
injection), mixing and combustion conditions at 
q=l. The Reynolds number based on the 
distance to the injectors for the flow compressed 

by the shock was - 0.5 *10 6 , which is consistent 
with that of a laminar boundary layer. The 
heat flux in this region agrees very well with 
laminar predictions. For the tare run, along the 
centerline (Fig 10a), the heat flux rises greatly 
downstream of the injectors, indicating tripping 
of the boundary layer by the injector holes. 
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(Without tripping, the boundary layer should 
remain laminar all along a model of this size.) 

Midway between the injectors (Fig 10b), the 
heat flux continues to drop downstream of the 
injectors until x = 32 cm, where it begins to rise 
somewhat. It is likely that this far 
downstream, the effect of the tripping of the 
boundary layer has spread in the spanwise 
direction sufficiently to effect the outboard 
gages. With injection, just downstream of the 
injectors (x = 23.5 cm) high heat flux levels are 
seen along the centerline, but not midway 
between the injectors. It is therefore likely 
that the plume shock does not reach the latter 
gage. Farther downstream (x = 26.5-35 cm), 
midway between injectors, large increases in 
heat flux occur upon injection. It is assumed 
that these areas must be in the compressed and 
heated flow aft of the plume shock system. 
Between x=29 cm and x=35 cm along the 
centerline, strong hydrogen film cooling is seen 
along the centerline, but not midway between 
the injectors. 

Surface Static Pressure 

Figure 11 shows the surface static pressure 
profiles for tare (no injection), mixing and 

combustion runs at q=l and combustion runs at q 
=2. Fig. 11a is for the centerline transducers, in 
line with an injector, while Fig. lib is for the 
two outboard rows of transducers, each midway 
between injectors. 

For the run with no injection, (circle symbols), 
all surface pressures were within a range of 
10%. Upstream of the injectors, the pressures 
for all runs were constant within 9%, showing 
no effect of the hydrogen plumes or the plume 
shock system. For the mixing and combustion 
runs at q = 1, substantial pressure rises were seen 
at the transducers 1.4 and 4.3 cm downstream of 
the injectors, showing the effect of the plume 
shock system. For these cases, the pressures at 
the transducers 10 and 12.9 cm downstream of 
the injectors had relaxed to the tare values. 
Also plotted in Fig. 11 are the CFD results for 
the combustion case. On the centerline, the 
calculated static pressure agrees with the 
experiment reasonably well. However, 
midway between the injectors the calculated 
values are slightly higher than those from the 


experiment. For the combustion run at q— 2, the 
pressure rises 1.4 and 4.3 cm downstream of the 

injectors were roughly twice those for q=l, and a 
noticeable pressure rise persisted to the center 
transducer 12.9 cm downstream of the injectors. 

From the wedge surface pressures alone, the 
difference between the mixing and combustion 

runs at q=l was quite small. Figure 12 shows a 
computation of the wall centerline static 
pressure for a mixing and a combustion case. 
The results show that for these conditions, the 
difference in pressure rise due to combustion is 
relatively small. This is due to the unconfined 
nature of the flow and the fact that the 
pressure upstream of the injector is low. Future 
tests in the facility will be operated at higher 
stagnation pressures, which will increase the 
pressure rise due to combustion. 

CONCLUSIONS AND FUTURE PLANS 

The results of the first propulsion related entry 
in the Ames 16 Inch Shock Tunnel were 
presented. Issues with regard to flow quality, 
flow contamination, test time and 
repeatability demonstrated the facility 
capability for high speed propulsion testing. 
Limited data on 30° flush wall hydrogen fuel 
injectors was obtained and the results were 
compared with a 3-dimensional Navier Stokes 
code. 

Future plans are to extend the operational 
pressure in the driver tube up to 680 Atm which 
will greatly enhance the facility simulation 
capability. Combustor inlet static pressure 
above 05 atm at conditions representing flight 
at Mach 12-16 will be obtained with this 
improvement. A new test cabin has been built 
which will be used to support testing of large 
scale models for NASP flowpath and NASA 
hypersonics research. 
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Table 1. Test Conditions 

Stagnation Pressure 
Stagnation Enthalpy 
Free Stream Velocity 
Free Stream Static Pressure 
Free Stream Temperature 
Free Stream Mach Number 

Fuel Total Pressure (q=l) 
Fuel Total Temperature 
Fuel Mach Number 


240 Atm 
95 MJ/kg 
3919 m/s 
.065 Atm 
1207 K 

5.8 

8.8 Atm 

300 K 
1 


Table 2. Seven-Reaction Model Arrhenius Rate Coefficients 


No. Reaction 


A 


N 


E, kcal/g-mole 


1 H2 + 02 <=> OH + OH 0.17E+14 

2 - H + 02 <=> OH + O 0.142E+15 

3 OH + H2 <=> H20 + H 0.316E+08 

4 O + H2 <=> OH + H 0.207E+15 

5 OH + OH <=> H20 + O 0.55E+14 

6 H + OH <=> H20 + M 0.221E+23 

7 H + H <=> H2 + M 0.653E+18 


0.0 

48.15 

0.0 

16.4 

1.8 

3.03 

0.0 

13.75 

0.0 

7.0 

-2.0 

0.0 

-1.0 

0.0 
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Figure 1. Ames 16-Inch Shock Tunnel 
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Figure 2. Test Hardware 


10 





Figure 3. Holographic Interferometry Apparatus 



Figure 4. Computational Approach 
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Figure 7. Interferogram for Case with Injection of Hydrogen into Air 



Figure 8. Shadowgraph of Case with Injection of Hydrogen into Air Compared with CFD 
Laterally Averaged Density Contours 
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ABSTRACT 

Flow characteristics of the 16-Inch Shock 
Tunnel at NASA Ames Research Center have 
been determined for purposes of providing 
hypersonic propulsion simulation capability. 
The key tunnel operating parameters are the 
incident shock speed and reservoir pressure and 
enthalpy. Flow characteristics of concern are the 
nozzle exit pressure, temperature, Mach number, 
Reynolds number, chemical composition, and flow 
uniformity. Surface mounted gages (for pressure 
and heat transfer) and nonintrusive optical flow 
diagnostics (emission and absorption spectroscopy 
and holographic interferometry) are used to 
verify tunnel conditions. Experimental 
measurements are used to validate computational 
anafysis for predicting facility performance, and 
CFD is used to interpret the free stream optical 
diagnostic measurements. 


* Research Scientist, Member ALAA 

** Research Scientist, Associate Fellow AIAA 
This paper is declared a work of the U. S. 
Government and is not subject to copyright 
protection in the United States. 


INTRODUCTION 

The 16-Inch Shock Tunnel at the NASA 
Ames Research Center has recently been 
identified as a viable facility for performing 
hypersonic propulsion research. It was designed 
and constructed in the early 1960's as part of a 
counter flow free-flight facility for purposes of 
simulating flight conditions for study of reentry 
aerothermodynamics. The facility was 
reactivated and upgraded under the National 
Aerospace Plane Technology Maturation Program 
to provide a test bed for large scale X-30 
flowpath testing at Mach numbers in excess of 12. 
This requires high dynamic pressure, high 
enthalpy flow corresponding directly to an 
airbreathing ascent vehicle trajectory. All Mach 
number simulation capabilities presently 
developed in the 16-Inch Shock Tunnel are with 
specific regard to an airbreathing ascent vehicle 
with a constant dynamic pressure of 1,000 psf. 

The operational characteristics of the 
shock tunnel are primarily determined by wall 
surface measurements of pressure and heat 
transfer. In the shock tube portion of the facility, 
the critical operational parameters are the 
incident shock speed, reservoir pressure, and 
enthalpy. These parameters are either measured 
directly by surface gages or computed using shock 
tube codes. Accurate determination is essential to 
ascertain flow characteristics, since reservoir 
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parameters serve as the initial conditions for 
computations. At the nozzle exit, flow 
characteristics of concern are pressure, 
temperature, Mach number, chemical 
composition, and flow uniformity. Wall surface 
mounted gages for pressure and heat transfer and 
optical diagnostics (emission and absorption 
spectroscopy and interferometry) are used to 
verify these tunnel operation conditions. 
Experimental measurements are used to validate 
computational codes for predicting facility 
performance. Conversely, instrumentation is not 
readily available for measuring all flow 
parameter profiles, so computational fluid 
dynamics must be used to compute the free stream 
properties. Similarly, CFD is used to interpret 
the * optical diagnostic measurements and 
calibrate some in situ probes. 

This paper will document the current 
capability of the 16-Inch Shock Tunnel by: 
confirming surface and flow field measurements 
with CFD computations; extending the 
understanding of the physics of facility 
operations by using computational models which 
allow the determination of system parameters 
not presently measurable; and performing data 
interpretation of flow field optical diagnostic 
results. Formal calibration of the shock tunnel 
flow properties with due regard to data precision 
and accuracy resulting from a thorough error 
analysis has not yet been completed. 

HYPERSONIC PROPULSION SIMULATION 
REQUIREMENTS AND CAPABILITY 

The use of continuous flow facilities for 
high speed propulsion testing is currently limited 
to flight Mach numbers of 12 or less. Simulation 
at higher Mach numbers requires the use of a 
pulse facility^ since it is only in these 
facilities that the high pressures, high 
enthalpies and relatively uncontaminated flow 
can be produced. A major concern in using pulse 
facilities for propulsion testing is the validity of 
the combustion data recorded during the 
relatively short test time. The test time in such 
facilities is on the order of milliseconds. 
Facilities of this type currently available are 
reflected shock tunnels and expansion tubes. At 
the higher Mach numbers, expansion tubes 3 ' 4 
may have an advantage over reflected shock 
tunnels in that that they may not stagnate the 


test gas, thus avoiding the production of a 
dissociated test flow. However, shock tunnels 
tend to produce longer run times, and currently 
operate at higher pressures. There are several 
types of shock tunnels, differing mainly in the 
method of operation of the driver tube. These 
types are the resistance heated driver 3 , the 
combustion heated driver^ 3 , and the free piston 
driver^. 

Proper flight simulation requires that 
the facility provide a combustor inlet flow with 
the suitable pressure, temperature, Mach number, 
Reynolds number, and test gas composition. Also 
of concern is that the test time be long enough for 
flow to become fully established^, which is on 
the order of the time required for several model 
lengths of test gas to pass through the combustor. 
Due to its method of operation and unparalleled 
large scale, the Ames 16 Inch Shock Tunnel has a 
unique capability for providing proper inlet 
conditions for high enthalpy (Mach 12-16) 
combustor testing using large scale test components 
combined with relatively long run time. 


FACILITY OVERVIEW AND OPERATIONAL 
PROCEDURE 

A schematic of the 16-Inch Shock Tunnel 
is shown in Fig. 1. Details of the driver tube 
configuration and operation can be found in Ref. 8. 
The driver section consists of a tube 70 ft (21 m) 
long with an inside diameter of 17 in (43 cm). The 
driven section is 85 ft (26 m) long with an inside 
diameter of 12 in (30 cm). The shock tunnel 
received its name for the 16-inch naval cannons 
used to construct its driver section. The shock 
tunnel is rated at 10,000 psia maximum driver 
pressure. The primary diaphragm is a flat, 304 
stainless steel plate, pre-scored to a depth of 
approximately 15% of its thickness. The 
diaphragm thickness is 1/16 in per 1,000 psia of 
driver pressure. The facility operates as a 
reflected shock tunnel with a thin sheet of Mylar 
separating the driven tube from the nozzle and 
test section. The driver tube is instrumented with 
pressure transducers at three axial locations 
evenly distributed along the length of the tube. 
The driven tube is instrumented with pressure 
transducers and shock detectors at five axial 
locations. 

The contoured Mach 7 facility nozzle is 19 
ft (5.8 m) long and has an exit diameter of 39 in 
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(0.99 m). Interchangeable throat sections are used 
to vary the nozzle area ratio. Recent results have 
been obtained for area ratios of 190, and a 
minimum of 95 is attainable without sacrificing 
test time or ideal shock tube end wall behavior. 
The nozzle contains ports at six axial locations for 
surface mounted instrumentation and optical 
access. The test cabin is a 72 in long by 54 in 
square cross section box located immediately 
downstream of the nozzle exit. It is designed 
with large doors on four sides providing easy 
personnel and optical access to test articles. 
Nozzle exit conditions are measured in this test 
cabin using a 35 probe pitot rake. Probes are 
located at one inch intervals and provide 
measurement of impact and static pressure and 
stagnation point heat flux. 

Operating conditions are achieved by 
combustion heating helium with a nearly 
stoichiometric mixture of hydrogen and oxygen. 
The gases are loaded into the driver tube through 
a 2 in inside diameter manifold which rests along 
the bottom of, and spans the entire length of, the 
driver tube. The manifold has small (0.020 in 
diameter) orifices evenly spaced along its length 
and is designed to promote uniform gas mixing. 
The loading procedure introduces the hydrogen 
first followed by the premixed combination of 
helium and oxygen. This combustible charge is 
ignited by four tungsten wires strung down the 
center of the driver and heated to approximately 
2,000 K by the discharge from a high voltage 
capacitor bank (C=150 |iF and V=14 KV). 
Combustion is complete approximately 25 ms 
after ignition, and the flat, prescored diaphragm 
breaks at the predetermined pressure. 

The shock speed, which determines the 
reservoir enthalpy and hence the flight 
equivalent Mach number, is varied by adjusting 
the ratio of the driver tube post-combustion 
pressure to the driven tube pre-test fill pressure. 
Thus, for the Mach 12, 14, and 16 conditions, the 
driver tube initial parameters remain identical. 
Driver tube operation has been recently 
demonstrated at pressures of 8,000 psia; however, 
efforts are being made to increase the operational 
pressure to the maximum rating of 10,000 psia. 

RESERVOIR CONDITIONS 

Reservoir conditions are those behind the 
reflected shock at the end of the driven tube. 
These conditions can be computed using standard 


shock tube relations. Variables in these relations 
are incident shock speed, initial fill pressure, 
measured reservoir pressure, and test gas 
composition. Equilibrium calculations are used to 
determine the conditions across the incident and 
reflected shocks. For a range of shock speeds, the 
shock tube can be operated in either the tailored 
interface mode or the non-tailored interface 
mode. For tailored interface operation, the 
reflected shock passes through the contact 
discontinuity without reflection. The reservoir 
conditions are those behind this single reflected 
shock. This is seen in reservoir pressure traces as 
a rapid rise to a constant pressure plateau. For 
non-tailored operation, shock or expansion waves 
successively reflect off the contact discontinuity 
and driven tube end wall causing a variation in 
the reservoir properties. These waves will decay 
quickly, usually after one or two reflections. This 
condition is marked by a more gradual settling to 
the constant pressure plateau. This method of 
operation, also referred to as the equilibrium 
interface 11 method, defines the beginning of the 
test time as that time after which reservoir 
properties cease to change significantly. For the 
equilibrium interface method of operation, if the 
calculated pressure behind the reflected shock is 
not equal to the measured pressure, an isentropic 
expansion to the measured pressure is performed 
to establish the equilibrium reservoir conditions. 
The test time is complete when the contact 
discontinuity reaches the reservoir and driver 
gases begin to enter the nozzle. Determination of 
this time period will be discussed in greater 
detail later in the paper. 

A graphical representation of this 
equilibrium interface generation process is shown 
in Fig. 2a (Ref. 21). The computed density 
contours are plotted on an X-T diagram along the 
full length of the shock tube. The computation 
begins with a simulation of the pressure 
distribution in the driver tube prior to diaphragm 
rupture and includes time dependent behavior, 
geometry variations, an approximation of viscous 
losses, and full chemical kinetics. The test gas 
region is shown at the far right of the diagram 
and in the expanded region of the plot. Evident is 
the incident shock generated upon primary 
diaphragm rupture and the several successive 
shock reflections establishing the pressure rise 
and the equilibrium interface condition. Figure 2b 
shows a comparison of experimental and 
computational driver gas pressure. The 
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computation which includes a sinusoidal pressure 
gradient in the driver tube agrees well with the 
measured pressure. Examining the results of this 
same computation for pressure at the driven tube 
reservoir also shows reasonable agreement with 
experiment (see Fig. 2c). This simulation 
capability allows for detailed examination of 
tunnel operating behaviors not previously 
available. Future expansion of the shock tunnel 
operational envelope will make use of these 
computations to optimize performance. 

Facility operation has been demonstrated 
at shock speeds ranging from 2. 6-3.5 km/s. 
Figures 3 and 4 show that these conditions match 
total enthalpies corresponding to Mach 12 to 16 
for the X-30 flight trajectory. Experimentally it 
is "found that, using nitrogen as the test gas, 
tailored interface operation is approximately 
achieved for a shock speed near 3 km/s (Mach 14 
enthalpy). However, preliminary data indicates 
that using air as the test gas, tailoring is not 
achieved for shock speeds over our demonstrated 
operating range. It is found that for these 
conditions the equilibrium interface is 
established about 2 ms after arrival of the 
incident shock. Reservoir pressures are nearly 
constant from this point until the arrival of the 
expansion wave, approximately 17 ms after 
incident shock arrival for present conditions. 


NOZZLE EXIT CONDITIONS 

The facility nozzle exit conditions were 
computed using the nonequilibrium nozzle code 
NENZF^. This is a standard inviscid 
nonequilibrium 1-D engineering code which has 
been updated with NASP standard chemical 
kinetic rate coefficients. Inputs to the code are 
reservoir conditions and nozzle contour. An 
effective nozzle area ratio of 170 is used in the 
sample calculations instead of the actual area 
ratio to account for viscous displacement effects. 
This was determined from experimental 
measurements of reservoir, pitot, and static 
pressures. Based on recent experimental results 
discussed below, the code appears to be 
satisfactorily predicting nozzle exit conditions. 

A recent experimental study of injection 
and mixing in hypersonic flows 13 provided an 
opportunity to measure free stream static 
pressure. The tunnel operating conditions used in 
this study are included in Table 1. For this test, a 


flat plate was installed at an 11 degree angle of 
attack to the flow and hydrogen was injected 
from the surface at an angle of 30 degrees. The 
model was inclined to the flow in order to 
increase the static pressure. For the highest 
tested driver operating pressure, 408 atm (6,000 
psia), the resulting reservoir pressure was 347 
atm (5,100 psia). At this reservoir pressure, it 
was seen that the static pressure as measured on 
the surface of the flat plate was just over 1/3 atm 
(5 psia). The measured pressure agreed with the 
pressure as predicted by NENZF to within 5%. 

More detailed computations using an 
axisymetric viscous code with finite rate 
chemistry and a quasi-one-dimensional code are 
complete and have been reported in Ref. 14. 
Comparisons of results from these two codes with 
NENZF results are shown in Figs. 5a, b, and c. 
Detailed computations predict gas temperatures 
at nozzle exit to be almost 30% lower than those 
computed using NENZF, while comparisons of 
static pressure and axial velocity are good to 
within about 6%. A lower O-atom recombination 
rate chosen for use in the detailed computations 
can account for the difference in temperature. The 
axisymetric flow code predicts a pressure 
gradient of about 6% across the inviscid core. A 
static pressure profile to be measured later in the 
program will provide a valuable validation for 
these computations. 

In order to achieve higher static 
pressures recommended for future combustor 
testing, the driver tube will be operated at 
pressures of at least 544 atm (8,000 psia). Table 2 
contains computed nozzle exit conditions for future 
tests to be run at Mach 12, 14, and 16 equivalent 
enthalpy at the 544 atm reservoir pressure 
condition. 

Another important nozzle exit parameter 
is the extent of the inviscid core flow. A previous 
calibration of the facility 6 measured the core 
flow at a position well downstream of the nozzle 
exit to be roughly 50% of the tunnel diameter. 
Correlations for boundary layer development 
indicate that for the new test cabin the core flow 
should be 70% of the diameter at the nozzle exit. 
Recent measurements with the instrumented pitot 
rake in the test cabin show a core flow of 
approximately 70% of the nozzle exit as 
predicted from earlier measurements (see Fig. 6). 
Uniformity of the inviscid core flow is also 
important. The pitot pressure profile shows a 
standard deviation in impact pressure of less 
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than 5%. As further evidence of uniformity. Fig. 

7 shows a shadowgraph of flow over individual 
pitot heads on the facility rake for a Mach 16, 
6000 psia driver pressure test condition. The 
shock structure on the pitot heads is clear, and 
there is no evidence of shocks from the nozzle 
disturbing the core flow. Plans for the calibration 
of the facility with the new test cabin and rake 
will be discussed in the conclusion. 

Turbulence simulation requirements for 
scramjet propulsion testing are not fully 
understood as of the present, and modelling of 
turbulent flows is not yet a mature technology. It 
is clear, however, that a fully turbulent boundary 
layer must be ingested by the combustor to 
simulate X-30 conditions. The free stream unit 
Reynolds number provided at nozzle exit for the 
8,000 psia driver pressure condition is determined 
to vary from 6xl0^/m to 2x10^/ m for the range of 
Mach numbers from 12 to 16 (see Table 2). Using a 
boundary layer transition criteria developed at 
Calspan^O anc [ found to be valid in their shock 
tunnel, the predicted transition point Reynolds 
number is about 1.5x10^. Given the large scale of 
the 16-Inch Shock Tunnel, a long enough inlet can 
be accommodated to provide ingestion of a fully 
turbulent boundary into the combustor. 

TEST TIME CONSIDERATIONS 

Though reservoir pressure may be nearly 
constant for a relatively long period, the 
composition of the gas passing through the nozzle 
is contaminated by driver gas much earlier. The 
arrival of driver gas contamination is dependent 
on several mechanisms including shock wave 
boundary layer interactions, contact discontinuity 
instability, and simple drainage of the test gas 
through the nozzle. With the current throat size, 
all of the test gas should drain from the reservoir 
in about 12 ms if there were no mixing of the 
driver and driven gases. This represents a 30% 
reduction in test time based on reservoir pressure 
alone. Mixing of the driver and driven gases at 
the contact discontinuity reduces the 
uncontaminated test flow even further. Several 
measurements have been used to determine the 
time of arrival of the helium (the primary 
constituent of the driver gas) including nozzle 
wall static pressure, heat transfer, free stream 
static pressure, and spectroscopic methods such as 
total radiation and laser absorption 


measurements^'^*. Similar conclusions are 
drawn from all of these measurements: initially 
there is a period of 2-3 ms dedicated to nozzle 
start up and decay of equilibrium interface 
transients followed by a nearly constant property 
flow period of 3-5 ms. As helium infiltrates the 
test gas, the nozzle wall static pressure and heat 
flux decreases (Fig. 8a and b). This is due to the 
increase in value of the ratio of specific heats of 
helium. Figure 8c shows the free stream static 
pressure as measured by a special static pressure 
probe adapter mounted on the pitot rake. This 
indicates that the nozzle start-up process is 
expanded in time as measured in the test section 
far downstream of the nozzle. Since helium and 
the other driver gas constituents are necessarily 
colder than the shock heated driven gas, the 
value of total radiation will also decrease as the 
helium contamination increases (Fig. 8d). The 
free stream OH temperature measured with a 
scanning laser absorption system^ is shown in 
Fig. 8e. (These temperature measurements are 
limited to at a few, low area ratio, nozzle 
locations.) For a few selected tests, the driven 
tube air was saturated with water. The shock 
heating process formed OH which was detected 
by the laser absorption diagnostic when 
expanded in the nozzle. Since detectable OH 
could only be present in the driven gas, the test 
time was assumed to be complete when the OH 
signal disappeared. The character and duration 
of the nozzle start up transient is well captured in 
the temperature measurements, and the constant 
temperture region is consistent with the test time. 
The unreduced OH absorption trace (not shown) 
illustrates the decay in driven gas OH mole 
fraction commensurate with the end of test time. 
An axisymetric, nonequilibrium nozzle expansion 
flow code*- 7 was used to compute the temperature 
at the first nozzle port. Fig. 8e shows the 
average measured temperature to agree with the 
computed temperature to within 1 percent. 

Nozzle wall static pressure is currently 
used on each run as the primary measure of test 
time. Recent results of these nozzle wall static 
pressure measurements are shown in Fig. 9. An 
initial transient period is indicated, followed by 
a period of approximately constant pressure, 
defined as the data averaging period. The 
indicated drop in static pressure correlates with 
the arrival of driver gas contamination confirmed 
by other methods of measurement, and is here 
used to indicate the end of the data averaging 
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period. Following this decrease in static pressure 
there is a period of nearly constant pressure, 
which is interpreted as flow of entirely driver 
gas. Therefore, three distinct test time periods 
can be defined: 1) a period of about 17 ms over 
which the reservoir pressure is nearly constant, 2) 
a period of uncontaminated test gas flow which 
includes the 2-3 ms for nozzle start up and the 
time to reach equilibrium interface, and 3) an 
averaging period of 3-5 ms in which the flow 
properties are essentially constant. 

Figure 10 plots test time as per definitions 
2) and 3) for different flight Mach numbers, and 
Fig. 11 plots test gas slug length. Slug length is 
defined by the expression, tU/L, where t is the 
facility test time, U is the flow velocity, and L is 
the characteristic length of the test article. 
These figures emphasize an important 
characteristic of the 16-Inch Shock Tunnel, 
namely its large temporal scale. When 
performing propulsion testing in pulse facilities, 
it is generally considered necessary that three 
gas exchanges through the combustor take place 
to assure flow establishment and hence proper 
assessment of mixing and combustion efficiency 
For facilities where test time is modest, the 
combustor model must be of a small scale to 
accommodate this criterion. For the Ames 16- 
Inch Shock Tunnel, however, with its long test 
times and commensurately long slug lengths (in 
excess of 10 m), virtually full scale X-30 combustor 
modules can be tested with confidence that the 
flow establishment criteria is sufficiently 
satisfied. 


CHEMICAL STATE OF TEST GAS 

As indicated in Tables 1 and 2, the 
reservoir gas temperature is between 4,600 and 
7,000 K. These high temperatures are 
characteristic of all reflected shock tunnels, and 
indicate a highly dissociated reservoir gas. The 
nozzle expansion process is nonisentropic and the 
test gas becomes chemically frozen^®. The result 
is a test gas which contains NO and dissociated 
oxygen. These species are known to affect 
combustion due to their influence on the ignition 
delay process and heat release^. For current 
operating conditions in Tables 1 and 2, the amount 
of NO in the test gas is about 5%. This is, 
however, not necessarily considered important in 
the scramjet combustion environment, since the 


reduction in ignition delay caused by NO will not 
significantly effect the mixing limited 
combustion process controlling propulsion at Mach 
number greater than 14. The dissociated oxygen 
entering a combustor, however, will produce a 
higher pressure rise due to combustion, since the 
recombination of atomic oxygen creates a mixture 
with an effective higher heating value. The 
amount of free stream oxygen depends, in general, 
on pressure, nozzle expansion rate, and nozzle 
length. Higher pressure promotes more collisions 
and hence higher recombination rates reducing O- 
atom mole fraction. More rapid expansion in the 
nozzle reduces collisions and allows higher mole 
fractions of O-atoms prior to the chemical 
freezing point. Longer nozzles produce less 
dissociated oxygen due to the longer residence 
time in the nozzle available for recombination 
collisions. The scale of the 16-Inch Shock Tunnel 
allows for a long nozzle (5.8 m) and a small area 
ratio (typically 170 for results presented here and 
as low as 95), resulting in relatively low mole 
fractions of NO and dissociated oxygen. Assumed 
is a constant reservoir pressure of 375 atm and 
values of A/ A* as shown specific to each Mach 
number. A plot of the variation of both NO and 
O-atom mole fraction versus Mach number and 
reservoir temperature is shown in Fig. 12. 
Although the NO mole fraction remains 
virtually constant for the considered range of 
Mach numbers at a level of about 5%, dissociated 
oxygen varies by more than an order of 
magnitude, from over 5% near Mach 16 to less 
than 0.2% at Mach 12. 


CONCLUSIONS AND FUTURE PLANS 

The NASA Ames 16-Inch Shock Tunnel 
has been successfully reactivated and its nozzle 
exit flow parameters characterized for purposes 
of hypersonic propulsion testing in the Mach 
number range of 12 to 16. An extensive 
experimental and computational effort have 
served to verify its applicability for X-oO 
flowpath testing as defined by the NASP 
program office. Nozzle exit static pressures are 
sufficiently high to provide combustor inlet plane 
pressures in excess of 0.5 atm (i.e. pressure 
achieved after the flow is processed by a cowl 
shock). The measured test time (3 to 5 ms) and 
test gas slug length (more than 10 m) coupled 
with its large nozzle diameter and core flow 
fraction allow proper flow establishment in full 
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scale combustor modules. Free stream Reynolds 
numbers and dissociated oxygen and NO mole 
fractions have been identified. Computational 
tools for shock tube analysis and nozzle exit flow 
quality have been developed and are serving to 
allow enhanced understanding of facility 
operation. 

An extensive calibration of the facility 
with its new test cabin will soon be completed. 
This calibration will include three different sets 
of conditions: Mach 12, 14, and 16 flows using air 
as test the gas. The primary purpose of this 
calibration process is to provide a matrix of 
documented conditions for X-30 combustor 
flowpath testing. The predicted conditions are 
given in Table 2. Throat sizes (hence, the nozzle 
area ratio) will, in general, be different for each 
condition. 

The primary goal of this calibration 
period will be to obtain experimental 
measurements including impact pressure, static 
pressure, and stagnation heat transfer at the 
nozzle exit rake. These results will be correlated 
with tunnel operating conditions computed using 
the measured quantities of shock speed and 
reservoir pressure to verify total enthalpy and, 
therefore, flight equivalent Mach number. 
Impact pressure measurements will be used to 
determine the extent of the inviscid core flow and 
the uniformity of the gas passing through the 
combustor. The static pressure probe data will be 
used primarily for determining test gas mass flow 
correlations, since this parameter is more 
sensitive than impact pressure to changes in 
reservoir conditions and in effective area ratio. 
The static pressure is also a sensitive indicator of 
test time. Included in this primary goal is the 
requirement of validating the computational 
fluid dynamics models that are being used to 
predict nozzle exit flows. Validation of these 
codes is essential to the future use of the facility 
given that all free stream parameters of interest 
cannot be measured experimentally. In fact, the 
static pressure probes used on the pitot rake 
require a computational calibration to provide for 
their quantitative use. This computational 
calibration effort is currently underway. 

All of these data will be critical for 
determining combustor inlet mass capture, a 
parameter essential for proper use of this facility 
as a propulsion test bed. The combustor planned 
for the first entry in the shock tunnel will be 
instrumented with inlet and cowl leading edge 
pressure transducers to assess the free stream flow 


parameters; however, CFD analysis will be 
required to mate the measured tunnel operating 
parameters with combustor inlet geometries and 
provide a more precise value of combustor inlet 
mass capture. 

It is important to note that although the 
present facility flow characterization and 
calibration effort is directed toward propulsion 
testing and research studies, the 16-Inch Shock 
Tunnel is not restricted to this use. Its will be 
valuable to experimental and computational 
research involving real-gas blunt body 
aerothermodynamics. * This includes flight 
trajectories for spacecraft that will be studied as 
part of the Mars mission program and NASA’s 
efforts to return to the lunar surface. Future plans 
will include calibration of test conditions 
required for these and other flight programs. 
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Flight Mach Number 


Mach 14 Mach 16 


Pressure, Atm 

245 

347 

Enthalpy, MJ/kg 

9.6 

12.2 

Temperature, K 

5750 

6850 


Pressure, Atm 

0.074 

0.12 

Temperature, K 

1280 

1720 

Mach Number 

5.7 

5.6 

Re/ m, 1CP*6 

1.5 

1.7 


Pressure, Atm 

0.27 

0.35 

Temperature, K 

1865 

2400 i 

Mach Number 

4.5 

4.4 

Re/m, 1(T6 

3.0 

3.4 


Table 1. Hypersonic Mixing and Injection Test conditions. 


Flight Mach Number 


Mach 12 Mach 14 Mach 16 


rveservuir v_uiiumuua 

Pressure, Atm 

517 

490 

463 

Enthalpy, MJ/kg 

6.6 

9.2 

12.2 

Temperature. K 

4650 

5730 

6930 

1 . *i ; 

Nozzle Exit Conditions 


Pressure, Atm 

0.16 

0.15 

0.13 

Temperature, K 

832 

1235 

1676 

Density, kg/m**3 

0.067 

0.042 

0.027 

Velocity, m/s 

3345 

3868 

4366 

Mach Number 

5.9 

5.7 

5.6 

Re/m. 10**6 

6.0 

3.4 

2.1 


Table 2. Conditions for 544 atm (8,000 psia) driver pressure operation. 



Figure 1. Ames 16-Inch Shock Tunnel 
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Figure 3. Total Enthalpy vs. Simulated 
Flight Mach Number 
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Figure 4. Incident Shock Speed vs. Simulated 
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Fig. 5. Nozzle flow code comparisons: 
Nozzle exit profiles 
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Fig. 9. Expanded and annotated nozzle pressure trace 


• Data averaging period 

♦ Data averaging period + transient 


Flight Mach Number 

Fig. 10. Test Time versus Simulated Flight Mach Number 
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Fig. 11. Slug Length of Test Gas versus Simulated Right Mach Number 


Mole Fraction of NO and O 
at Nozzle Exit vs. Stagnation Temperature 


O 

- - NO 


Po=375 atm 
Nozzle Area Ratio=170 


stagnation 


Fig. 12. NENZF computations of NO and O-atom mole fractions versus 
stagnation temperature and simulated flight Mach number. 
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Abstract 

Current and future calculations of nonequilibrium 
shock layers require the use of a very large number of 
equations, due to a multiplicity of chemical species, 
excited states, and internal energy modes. The com- 
putational cost associated with the use of standard 
implicit methods becomes prohibitive; it is, there- 
fore, desirable to examine the potential of several 
methods and determine if any can be projected to 
be more efficient and accurate for large systems of 
equations. In this paper we examine the performance 
of several implicit schemes on some simple practical 
examples of reacting flows. The Euler equations are 
solved by three different implicit methods, and two 
methods of coupling between the fluid dynamics and 
the chemistry are studied. Several cases of stiffness 
are considered, and both one and two-dimensional 
examples are computed. We conclude on with some 
remarks on the accuracy, stability and efficiency of 
these various methods. 

I. Introduction 

The modern Computational Fluid Dynamics 
(CFD) tools are becoming increasingly useful in com- 
puting complex flow conditions, which generally in- 
clude non-equilibrium phenomena. There is a gen- 
eral need for increasingly complex modeling of the 
thermo-chemical properties of the gas, and for the 
modeling of larger systems. For example, the model- 
ing of shock layers around ablating bodies requires a 
very large set of chemical species and chemical reac- 
tions. Although some approximate formulations can 
be used in the preliminary design phase of space ve- 
hicles or experiments, the modeling of the complete 
kinetics is desirable or even required when the non- 
equilibrium effects become dominant: this happens 
for example as the flow expands around the shoul- 
der of a vehicle, or when the object is reduced in 

’Senior Research Scientist, member AIAA. Mailing ad- 
dress: NASA Ames Research Center, MS 230-2, Moffett Field, 
California 94035. 
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size for insertion into an experimental facility. The 
situation can be further complicated due to the fact 
that most flows realized in ground-based experimen- 
tal facilities are themselves not in chemical or ther- 
mal equilibrium. Another example concerns highly 
ionized and radiating flows, which are likely to be 
found at high re-entry velocities, or their equivalent. 
It may be required, for these cases of plasma condi- 
tions, to account for non-Boltzmann distribution of 
the excited states. This problem may also require 
us to model the plasma with a complete collisional- 
radiative model of the plasma, and convect all the 
excited states, forcing us to use a large number of 
equivalent species. In addition, the internal relax- 
ation processes will be locally stiff and difficult to 
model. These upcoming challenges in CFD technol- 
ogy will require the development of efficient meth- 
ods for a very large number of species, and for pos- 
sibly stiff couplings to complex internal processes. 
Since we want a method that allows us to reach 
the steady state with minimal computational effort, 
it seems desirable to use an implicit method. On 
the other hand, since we may need to couple the 
hydrodynamics to several other physical processes 
(collisional-radiative processes, radiation transport, 
electro-magnetic couplings, etc..), we may favor the 
use of the Operator-Splitting (OS) method. The lat- 
ter must be contrasted with the Fully-Coupled (FC) 
approach, which attempts to provide a more accu- 
rate and more stable way to couple the different pro- 
cesses. It is not clear a priori which method is more 
accurate, stable, efficient, or practical; the use of one 
method versus the other may depend on the type of 
flows being computed, the type of computer archi- 
tecture used, or even the personal preferences of the 
modeler. During the course of the present work, we 
will compare the FC and OS methods (for chemistry) 
and therefore add fuel to the debate. Our search for 
an efficient numerical method, extendable to large 
systems, will also include the effect of required grid 
accuracy on the solution, and its impact on the effi- 
ciency of the numerical approaches used. 


II, Numerical Methods 

II- A: The N x N Block-Tridiagonal method 

The Euler equations describe the convective pro- 
cess, and are written (in ID) as: 



where E } H are respectively the total energy and to- 
tal enthalpy, per unit volume. The internal energy 
is Ei = f T dT'C v} (T'), and the gas mixture 

follows the ideal equation of state: 

P = NkT = (7 - 1 )Ei = (7 ~ 1 )(£ - \pu 2 ) 

The individual species densities are denoted by p 5) 
and p — Ps- This formulation is for a single fluid 
(one mass-averaged velocity), in thermal equilibrium 
(one temperature). This equation is discretized over 
a finite size mesh to yield the form: 



The subscripts z, i ± ^ indicate that the variables 
are evaluated at computational cells (center) i, and 
cell interfaces i ± \ . V and S are respectively the 
cell volumes and surfaces. This finite-volume for- 
mulation will be used throughout the paper. The 
equation has been discretized in time as well, and 
the expression AQ describes the difference between 
the flow variable evaluated at two time levels, n and 
n -j- 1. The expression on the RHS of equation (2) 
must be further specified: the fluxes are a function of 
Q, and can be evaluated at time level n *f 0 through 
the linearization approximation: 

+0A(Q)(AQ) i+ * (3) 

where A = is the Jacobian matrix. The explicit 
Euler system of equations is obtained for 9 = 0, the 


implicit system for 9 = 1, while second-order time 
accuracy is obtained for 9 = 1/2. Second-order spa- 
tial accuracy is achieved by evaluating the fluxes at 
the cell interfaces i ± 

F <±i = 2^ F * + Fi±1 ) 

A final modification to the fluxes is made to assure 
monotonicity. The Euler system is an hyperbolic sys- 
tem, and has a set of real eigenvalues (characteristic 
speeds). The Jacobian can then be written in the 
form: 

A = T" 1 A T (4) 

where the matrix of eigenvalues 



is diagonal and real only, T, T" 1 are transfer matri- 
ces between the space of primitive variables Q and 
‘characteristic' variables, and c is the speed of sound. 
The spectrum of eigenvalues can be split into and 
positive and negative values, which indicate the di- 
rection of flow of the characteristic variables at the 
cell interface. The flux at an interface can now be 
written as: 

F (n+9) ^ F (") i+ 0 A + (Q)AQ,+^A-(Q)AQ i+1 (6) 

* + T 3 

where 

A ± = x- 1 • A ± • T 

and A* is the set of eigenvalues which are respec- 
tively positive (negative), zero otherwise. Using 
this formulation, the discretized version of the Eu- 
ler equations becomes: 

[—OAtAf t ] AQ,*_ i 

* f 

+[1 + 9&t A+ , - 0A*A* , ]AQ, (7) 

*13 1 5 

+[«A1A- j]AQ, +I = -F« 

The RHS of equation (7) can be modifed for mono- 
tonicity, while conserving its second-order accuracy 
in space. The technique used throughout this work 
follows closely the TVD method of Harten [l] . 
Greater stability is generally obtained if the implicit 
LHS of equation (7) has its spatial accuracy reduced 
to first order. This consists in evaluating the Jaco- 
bian matrices at the cell centers, according to the 
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characteristic flow direction. The final version of the 
system is obtained for the fully implicit case (0 = 1 ): 

[-AtiA*]AQ,-_i 

+[1 + At At - A“ A<] AQ, ( 8 ) 

+[A- + 1 A*]AQ <+1 = RHS^ 

This is a tridiagonal system of N x N block ma- 
trices, where N = N s + D + 1 , N s is the number 
of species, D is the spatial dimensionality. Solv- 
ing this system can be done by the standard tech- 
nique of gaussian elimination and back substitution, 
with LU decomposition of the block matrices (see 
for example [ 2 ]). This requires that the matrices 
that compose the diagonal band be inverted twice 
for each grid point. It turns out that the cor- 
responding algorithm has a number of operations 
that scales as N 3 . This technique is called the 
N x N Block-Tridiagonal solver, and will be used in 
this paper. 

A similar approach can be used for 2-dimensional 
flow, leading to a pentadiagonal system of N x N 
bock matrices. Another approach, which consumes 
less memory and has a lower number of operations, 
uses the technique of dimensional splitting [3]: we 
effectively solve the tridiagonal system twice, once 
for each direction^. This is the approach used here. 

II-B: The Scalar Tridiagonal method 


It is clear that as the number of species grows, 
the iV 3 dependence of the computational load will 
rapidly make this technique impractical. It is de- 
sirable then to search for a method that scales less 
rapidly with the number of species. One such method 
can be easily obtained by considerable simplification 
of the matrix structure. Note that the split Jacobians 
are bounded respectively from below and above: 

A+ =T" 1 .A + .T<max{A}l (9) 

A - = T - 1 .A~.T > min{A} • 1 (10) 

where we have used the (signed) maximum and min- 
imum eigenvalue present in A. The Jacobian matri- 
ces, thus replaced into equation ( 8 ) are proportional 
to the unit matrix. There is only a scalar opera- 
tion to perform, instead of a full block matrix in- 
version. The number of operations now scales as N: 
this scheme is called the Scalar-Tridiagonal solver, 
and will be compared to the previous one. 


fThis splitting is a form of Operator-Splitting. 


There are several disadvantages to the scalar 
technique: the first, loss of time accuracy, is not of 
immediate relevance, since we are mostly concerned 
with the achievement of steady-state. The second 
is a loss of accuracy: this is specially of concern in 
subsonic regions, where the spectrum of eigenvalues 
is originally very different from the maximum (min- 
imum) value. In supersonic or hypersonic regions, 
this is not a problem, since u >> c, and A a til, i.e. 
the spectrum is nearly scalar. We may expect there- 
fore some loss of accuracy, or even stability, when 
using the scalar method. 

II-C: The N s -Split Tridiagonal method 


We will also investigate another method, based on 
the formulation of a multi-fluid system of equations. 
Let us consider the following system of equations: 



up s 

p, + p,u 2 

uH, 


= 0 ( 11 ) 


and similar systems for other species. In this for- 
mulation, each species is attributed its own momen- 
tum density and energy density. In the limit of very 
strong coupling between the momentum and energy 
densities of each fluid component, we can enforce 
a unique velocity and unique temperature for this 
multi-fluid description. If we were to solve each sub- 
system by the block tridiagonal method, we would 
require the inversion of a 3 x 3 (4x4 in 2 dimensions) 
matrix. We repeat the method for each species, and 
the overall cost now scales as N s . For a large num- 
ber of species, we expect a cross-over between this 
method and the N x N block method, by compar- 
ing the costs; for example, 4 3 A r a versus (N s + 3) 3 . 
However, a further reduction in cost can be achieved 
with the following approximation. Assuming that 
all species have nearly equal molecular masses, and 
that their individual specific heats are nearly equal, 
we can replace the partial pressure: 

p , = n,RT~ P —P 
P 

and use a constant average 7 in the formulation of 
the derivatives which compose the Jacobians. These 
Jacobian matrices become then identical. This has 
a rather drastic effect: the block matrix inversions 
need to be performed once, instead of once for each 
species. This lowers considerably the overall CPU 
requirement. This formulation of the solver will be 
called Nj-Split Tridiagonal solver, and will compared 
with the two previous ones. 
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II-D: Hydrodynamics-Chemistry Coupling 

In the equations considered so far, the source 
terms on the RHS are non existent; for a reacting 
gas, there will be a source term W which operates 
on the species densities only. The total energy now 
includes the energies of formation: 

E=E t + l -pu 2 + Y. e ° 

3 

and is not not affected by the chemical reactions: 
the internal energy is obtained from the conserved 
total energy, after subtracting the kinetic and for- 
mation energies. The change in formation energy of 
the mixture, due to chemical reactions, will thus be 
converted into a change of temperature. 

The chemical source term can be treated im- 
plicitely as well: if 17 is the corresponding Jacobian 
for the source term (17 = ^y), solving for the chem- 
istry alone would read as: 

[1 + fiiAijAchQ = w[ n) Ai (12) 

By solving separately for the convective and chemical 
terms, one obtains two variations at the end of the 
computational step: the global variation will then 
be a direct sum of the the contribution from each 
process. 

AQ = A cv Q + A ch Q (13) 

This procedure is called^ the Explicit-Coupling 
method (EC). 

We see that in effect, we have split the com- 
putation in two parts, for each physical process. 
For that reason, this method is also called the 
Operator-Splitting (OS) method (see for example 
[3]). Another form of operator splitting consists in 
using the change induced by one process as a starting 
point for the other process: a temporary solution Q 
is used, such that: 

Q = Q (n ) + A CV Q (14) 

Q^+D = Q + A ch Q (15) 

where now the change induced by chemistry is ob- 
tained by using the modified solution Q in the ex- 
pression of the source term and Jacobian, W,17 in 
eq. (12). This formulation of the Operator-Splitting 
method is based on fractional steps, and is best de- 
scribed in [4]. Since we will examine both meth- 
ods, we will reserve the term Explicit-Coupling (EC) 

1 Also called the Loosely-Coupled approach. 


for the method described in eq. (13), and the term 
Operator-Splitting (OS) when using the method de- 
scribed by eqs. (14) and (15). 

Another approach is to solve for both the con- 
vective and chemical processes simultaneously. The 
chemical Jacobians can be brought into the LHS, and 
equation (8) is modified to: 

[-A+iAfJAQi-i 

+ [1 + Af At - A“ A t + 17, A*]AQ t (16) 

+[A- +1 A<]AQ, +1 = RHS™ 

where now the RHS includes the evaluation of the 
chemical source terms at time level (n). This method 
is called the Implicit, or Fully-Coupled (FC) ap- 
proach. The O matrix is dense, and the FC approach 
described above is possible only when combined with 
the N x N Block- Tridiagonal solver. Including it in 
the iV.,- Split solver would require serious modifica- 
tions, and has not been attempted here. Similarly, 
by approximating the Jacobian 17 by a scalar (using 
again the maximum eigenvalue), one could use the 
FC approach with the Scalar Tridiagonal solver. We 
found that in many cases this approximation usually 
leads to very poor results for the chemistry, and will 
not be used here, 

II-E: Chemistry Sub-Cycling 

There are additional modifications one can make 
when using the OS or EC approaches: since the fluid 
dynamics and chemistry are computed separately for 
a global time step At , one has considerable flexi- 
bility in the methods used for each process. No- 
tably, the accuracy of the chemistry can potentially 
be improved by sub-iterating (more precisely sub- 
cycling) the chemistry by using smaller time steps 
8t. This may be required to improve the accu- 
racy, because the chemical reactions are non-linear 
processes: linearization errors become important in 
some highly non-equilibrium situations. The cou- 
pling of the chemistry to the temperature can also be 
estimated at each sub-step, by looking at the change 
induced in the average formation energy of the mix- 
ture. When the chemistry is sub-iterated (SI), the 
global variation is obtained by using eq. (13), but 
when the global change due to chemistry is obtained 
as follows, using sub-iterations (m = 1, 2, . . .) of the 
chemistry: 

A ch Q (m+1) = A ch Q (m) + [1 + W St (17) 

Finally, the coupling of the chemistry to the con- 
vection can also be computed at each sub-step. For 
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example, after computing the variation A CV Q using 
one of the Tridiagonal solver listed above, one can lin- 
earize it during the global time interval. The global 
variation (for both processes) is then obtained as fol- 
lows: 

AQ( m+1) = AQ fm) + [1 + ft6*r l W St (18) 


This formulation replaces equation (13). We will 
reserve the term Sub-Iterated Coupling (SIC) to de- 
scribe this particular form of the splitting method 
between chemistry and convection. 

II-F: Performance 

The methods used can be classified, according 
to the treatment of the convective process and the 
method of coupling with the chemistry: the designa- 
tion of the methods which will be studied are listed 
in Table 1. 

The relative performance of all schemes is demon- 
strated in Figure 1. Figure 1-a (top) shows the CPU 
spent (per iteration and per grid point) by an im- 
plicit method, normalized by the same quantity for 
the explicit method. Since the explicit method scales 
almost exactly as the number of equations, both the 
Scalar and the N s -Split method will show a nearly 
flat behavior when plotted versus the number of 
species. This is confirmed in Figure 1-a. Notice also 
that the relative cost of the Scalar method is very 
small, while the N x N Block method climbs very 
rapidly: the latter is still quite expensive, even for a 
small number of species. The leftmost data point at 
N s = 5, for example, shows that the N x N Block 
method is 10 times more expensive than the explicit 
method. Although this number is not an absolute, 
and can be reduced after a strong effort in code writ- 
ing (by ‘hard- wiring* the operations, for example). 
At best, this time may be reduced by a factor of two. 
Still, the conclusion is inevitable: as the number of 
species grows, the implicit scheme is efficient only if 
it can be operated at large CFL numbers. Practically 
speaking, stability limitations will limit the CFL to 
the neighborhood of 4-5. These stability problems 
arise from transient phenomena, dimensional split 
errors and/or coupling errors with the chemistry or 
other internal relaxation processes. Higher values of 
CFL number can potentially be achieved when the 
flow is very close to the steady state and when the 
flow is non-stiff. Since we are mostly interested into 
reaching the steady state (and having to go through 
the transients) and into stiff problems, this is of lit- 


tle interest to us. These limitations will be demon- 
strated on some practical sample cases in the next 
sections. 

ILL Qne-Pimonsional Shock 

As a first test case, we will model the propagation 
of a 1-dimensional shock, from an impulsive start. 
This case will mimic the establish ement of a two- 
dimensional shock layer around a blunt body. We 
use a grid of 200 cells, evenly spaced, with a per- 
fectly reflecting wall on the right hand side. The 
flow is incoming from the left at high velocity, and 
impinges on the wall. A shock is created at the re- 
flection and propagates back upstream into the hy- 
personic flow. Although strictly speaking this flow 
is unsteady, the profiles become steady in a frame 
attached to the shock. The gas is air, composed of 
5 species (AT, O, N 2 , 0 2 , NO), the free stream Mach 
number is = 15, the free-stream temperature is 
= 300° A*. Three cases of free-stream pressure 
will be considered, leading to three stiffness condi- 
tions: 

easel : Poo = 10“ 5 atm 

case 2 : Pco = 10“ 4 atm 

case3 : Poo = 10“ 2 atm 

The stiffness is defined as the ratio of the largest 
time scale (here presumably the convective one) to 
the smallest (chemistry). The convective time scale 
is obtained from the choice of Courant (CFL) number 
we choose to run the simulation at. The chemistry 
time scale can be defined in two ways: 

• an intrinsic time scale, obtained from the maxi- 
mum rate of change of any chemical specie. For 
example, the chemical time scale will be the time 
required for a specie molar fraction to change by 
more than 10%, provided it is not close to zero 
initially. 

• a coupling time scale, defined as the time re- 
quired for the chemistry to modify any flow vari- 
able by (say) more than 5% . Since the chemistry 
affects mostly the temperature, this is the vari- 
able used in that case. 

The second time scale provides a global limitation on 
the time step to be used: if t he chemical effects dra- 
matically change the formation energy of the mixture 
during the time step, and if this Se° is large com- 
pared to the internal energy, the numerical solution 
becomes rapidly unstable. This has a profound effect 
on the choice of numerical methods to be used, for 
example, in combustion. In the remainder of this 
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paper, we always limit the global time step such 
that the estimated relative change in temperature, 
induced by chemical reactions, is smaller or equal to 
5%. We allow the use of large (CFL> 1) global time 
steps provided this condition is satisfied. Let us em- 
phasize that this restriction still allows us to consider 
stiff problems, where the stiffness is defined by using 
the intrinsic chemical time scale. Chemical equilib- 
rium can be reached rapidly, without significantly 
modifying the temperature; the flow conditions sim- 
ply must be such that the equilibrium values are not 
very different from the initial values, or that the en- 
ergies of formation are relatively small compared to 
the internal energy. 

The results presented in this section are obtained 
using the following methodst : 

1 FC/Block tridiagonal 

2 OS/Block tridiagonal 

3 EC/Block tridiagonal 

4 EC/Scalar tridiagonal 

5 EC/N S — Split tridiagonal 

The chemistry is always computed with a single it- 
eration. 

Let us look first at case 1, for Pea — 10 -5 atm. 
The profiles of temperature are shown (left scale on 
the plots) as well as the mole fractions of N and O 
atoms (right scale on the plots). Figure 2 shows these 
profiles for a calculation at CFL=2 and 4. Figure 2-a 
shows only the results for the FC/Block, EC/Scalar 
and EC/N, -Split methods: the remaining cases of 
EC/Block and OS/Block would show profiles exactly 
identical to the FC/Block method. The agreement 
between the other methods is also quite good. This 
is also true for the CFL=4 case (Figure 2-b, bottom), 
although to a lesser extent: in this figure, the curves 
for the EC/Block case are omitted, since they are 
identical to the FC/Block results. It appears there- 
fore for this case that the EC and OS methods are 
as accurate as the FC method. The N, -Split method 
shows slight errors in species concentrations near the 
shock, in the region of highest concentration gradi- 
ents, which worsen as the CFL number grows. The 
Scalar method has an overshoot at CFL=2, and can- 
not be operated at larger CFL numbers. All methods 
fail for larger CFL values. 

Figure 3 shows the same profiles (obtained when 
the shock reaches the same position) for a slightly 

tThe notation used has been mentioned in the previous 
section, and is summarized in Table 1. 


stiffer problem (case 2). Again, we had perfect 
agreement between the OS/Block, EC/Block and 
FC/Block methods, and the EC/Block profiles were 
omitted for clarity. The overshoots in mole fraction 
near the shock, for the Nj-Split method, are worse 
for this stiffer case. Again, the Scalar method works 
reasonably well for CFL=2, but fails for larger val- 
ues. These calculations were done using the standard 
mirnnod limiter in the convective fluxes, as described 
by Harten [3], with an entropy fix € ~ 0.1. When us- 
ing a more compressive flux-limiter, such as the ‘Su- 
perbee’ limiter, the calculation could proceed as well, 
although with very slight oscillations. Reducing the 
entropy parameter to e ~ 0.01 would lead to more 
severe oscillations. Therefore the rapid elimination 
of the transients can be best achieved by ensuring 
that sufficient numerical diffusion is present. The fi- 
nal flow solution therefore would need to be further 
sharpened, when the steady state is nearly achieved. 

Figure 4 shows the stiffest case for CFL=2. All 
methods failed for larger CFL values. It is remark- 
able that the FC/Block method failed for this case, 
while the EC/Block method gives the best results. 
The OS/Block method (which uses the fractional 
step approach) gives very similar results, and can be 
considered as accurate. Surprisingly, the EC/Scalar 
method is stable, although the species profiles show 
an unphysical kink in the relaxation region. In order 
to better determine which method is more accurate, 
we computed the same case on a larger grid (2000 
points) using the FC/Block method. By increasing 
the grid density, we achieved a ten-fold reduction 
of the stiffness of the problem. The FC/Block was 
then run successfully, and gave a very short relax- 
ation zone (see Figure 4b). We also attempted to 
better reproduce this relaxation on the coarse (200 
points) grid by either 1) sub-cycling the chemistry 
or 2) reducing the time step. Figure 4b shows the 
comparison, for example between the FC/Block re- 
sults computed on the high-density grid, with the 
EC/Scalar with 10 sub-iterations of the chemistry 
and an explicit calculation (CFL=0.2). The two lat- 
ter cases are not very different from the results of 
Figure 4a, i.e., neither the sub-iterations nor the time 
step reduction greatly improve the solution. It seems 
that all methods tend to overestimate the length of 
the chemical relaxation zone in stiff cases, although 
the final equilibrium result is accurate. We must 
conclude also that the EC or OS methods are more 
stable than the FC method in stiff cases: we also ob- 
served this feature on other stiff cases. The mixing 
of non-diagonal elements in the global Jacobians, be- 
tween convective and chemical terms, may make the 
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matrices more prone to ill-conditionality, and reduce 
the stability. 

It seems therefore that only relatively small val- 
ues of the CFL number can be effectively used for the 
transients, and therefore only the Scalar method, so 
far, is efficient. However, it is not accurate enough 
when the chemistry is very stiff. It also appears that 
the A^-Split method, in its current form, suffers from 
unphysical numerical species diffusion in the region 
of strong gradients, and for large time steps (this er- 
ror is inexistent in the explicit regime). Since there 
are other cases where implicit methods can have a 
significant impact, we will look also at expanding 
flows in the next section. 

IV. One-Dimensional Nozzle 

We will model a converging/diverging nozzle, 
with 150 grid points in the axial direction: the calcu- 
lations were performed using two-dimensional codes, 
and the grid used 10 points in the radial direction. 
Since we were interested only into axial profiles, this 
was considered sufficient for our purposes. Notice 
that now there will be an additional error due to 
the dimensional split in the implicit methods. The 
left boundary condition and initial state considered 
a gas at a pressure of 4.205 atm and 1000°K. How- 
ever, the gas composition was arbitrarily set to non 
equilibrium values by increasing the amount of disso- 
ciation: this had the effect of stiffening the chemistry 
in the subsonic region of the nozzle. The calculations 
were always started impulsively, and run at CFL< 1 
(explicit) until the shock exits the nozzle, before the 
implicit models were tried. We used this case to eval- 
uate the effect of sub-iterations and sub-coupling in 
the chemistry. 

Figure 5 shows the comparison of residual his- 
tory for the 4 implicit methods used, i.e. FC/Block, 
EC/Block, EC/Scalar and EC/N 5 -Split, and with- 
out any sub-iterations or sub-coupling. The resid- 
ual of the subsonic zone (solid line) and supersonic 
zone (dotted line) have been shown separately. A 
first break point in the curves shows the end of the 
explicit pre-calculation, used for elimination of the 
shock from the nozzle. The implicit scheme is then 
used, with a constant CFL=1.5, until a time of 5 mil- 
liseconds. At that point, the calculation is pursued 
further for the supersonic region only, the subsonic 
region remaining frozen. The residual shown is for 
the total energy, and is averaged over the entire vol- 
ume of the region considered. 


Figure 5-a shows the results for the FC/Block 
method, Figure 5-b for the EC/Block method. Both 
show good convergence properties, with a slight im- 
provement for the EC/Block method. The V^-Split 
method (Figure 5-c) has even better convergence 
properties, but the Scalar method (Figure 5-d) shows 
a non- vanishing residual for the subsonic region. All 
methods converge rapidly in the supersonic region 
when computed separately, indicating that most of 
the problems (if any) are located in the subsonic re- 
gion. 

The solutions obtained at 5 milliseconds are plot- 
ted in Figure 6, for the atomic oxygen mole fraction. 
The solution for the scalar method is slightly in error 
in the subsonic region, but quickly recovers during 
the expansion and leads to the correct final value. 
The Nj-Split method has the opposite behavior, i.e. 
has an error increasing with the distance along the 
nozzle: it seems that the species convection suffers 
from some unphysical diffusion of species, also no- 
ticed in the results of the previous section: there 
is a phase error between each species convection, 
which is irreversible. By contrast, the Scalar method 
correctly propagates the species, but does not ac- 
curately couple the convection with the momentum 
and energy equations. This may lead to fluctuations 
in pressure or temperature, which quickly disappear 
when the flow becomes near supersonic. 

The use of sub-iterations in the chemistry did 
not change the results for this case. Increasing the 
stagnation pressure and the stiffness slowly lead to 
noticeable effects. The most dramatic differences 
between the cases of sub-iterated and non-iterated 
chemistry were observed for very stiff systems, at 
the limit of stability. In order to demonstrate the 
effect of sub-iterations, or sub-cycling of the chem- 
istry, we consider a high pressure (400 atmospheres) 
case, with an initial temperature of 6000 °K , and 
a highly non-equilibrium initial composition (non- 
dissociated air). A constant time step of 5 x 10” 8 
seconds was assumed. Figure 7 shows the results 
of the chemistry integration (no fluid dynamics) for 
both non-iterated and sub-iterated (20 cycles) cases. 
It is clear that a single step of the chemical integra- 
tion, with At = 5 x 10~ 8 sec, leads to very large 
changes in species concentrations and temperature. 
This will significantly affect the remainder of the his- 
tory of the chemical integration. If the time step is 
not too large, the correct equilibrium values may be 
obtained in the final steady state: if the time step 
is large enough, unphysical values (i.e. negative con- 
centrations) may be obtained during the first step, 
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and the correct solution cannot be recovered. By 
using a smaller time step (10“ 8 sec), both methods 
give essentially the same history, and the same fi- 
nal values, which agree with the values obtained in 
Figure 7 for the sub-iterated case. It is clear then 
that in some severe cases, the sub-iteration of the 
chemistry can yield a higher stability and a higher 
accuracy: these cases may be found for example in 
high-pressure shocks, detonations, or strong ionizing 
shocks, when the grid used is coarse. 

The method of sub-coupling was also tested on 
some other stiff cases. It was found that the stabil- 
ity was slightly reduced when the sub-coupling was 
incorporated. When the chemistry is sub-iterated 
and sub-coupled, the convection of species during the 
sub-step St is estimated and included in the varia- 
tion. While this process accounts for the effect of 
the convection on the chemistry, it fails to account 
for the reverse process, and it fails to account for 
the influence of other convective terms, specifically 
the pressure waves. The coupling of the chemistry is 
not ‘in phase 1 with all the convective equations. In 
subsonic regions and behind shock waves, the pres- 
sure waves are a dominant process, and a signifi- 
cant error is made. We therefore recommend that 
no sub-coupling (SIC) be used if the chemistry is 
sub-iterated. 

V: Two-Dimensional Shock Layer 

The final test will be done for a 2-dimensional, 
axi-symmetric flow around a blunt body. This is a 
typical flow configuration of interest. The flow is 
modeled with a 144x80 grid, the free stream Mach 
number is A/oo = 25, the free stream tempera- 
ture and pressure are Too = 241.75°K and P ^ = 
1.65 10 -4 atm; the free stream is air, modeled using 
11 species N t O,N 2 , 02 ,NOtN+, 0 + 9 N} t Oi,NO+ t e- 
and a 15 reaction set from Dunn & Kang [5]. The 
blunt body shape is taken from the Apollo space- 
craft. 

The calculations were proceeded with the 
FC/Block method and the OS/Scalar method. Since 
the chemical changes were quite important in the 
shock layer, our restriction in time step due to the 
temperature changes (relative change < 5%, see sec- 
tion III) prevented us to compute the flow implicitely. 
Any attempt to increase the time step, and therefore 
to allow a larger change in temperature due to chem- 
ical effects led very quickly to flow instabilities. The 
comparison presented below is therefore between a 
FC/Block implicit method run at small time steps, 


and an Operator-Split method where the fluid dy- 
namics are computed explicitely. The pressure be- 
hind the shock is close to 0.5 atmospheres, with a 
temperature between 8,000 and 12,000 °K, and the 
chemistry is rather stiff, especially due to the reac- 
tions involving electrons. Figure 8 shows the compar- 
ison between the two methods along the stagnation 
line, and the agreement is very satisfactory. 

There is a lack of resolution of the shock front, 
and we proceeded to improve the results by adapting 
the grid in the neighborhood of the shock. Several 
adaptions were performed, first on the temperature 
gradients, then on the chemical gradients (A^)- Af- 
ter each adaption, the flow was computed further 
until convergence. The adaption procedure used the 
SAGE code [6] developed at Ames, and affected grid 
points in a direction normal to the blunt surface only. 
Figure 9 shows the comparison between the origi- 
nal, non-adapted solution and the results from the 
final adaption. Since the results from the FC/Block 
and OS/Scalar methods were found to be agreement, 
and since the latter method is considerably more ef- 
ficient, only the OS/Scalar method was used for the 
adapted cases. We see in Figure 9 that the peak 
temperature has changed significantly (15 %), and 
so has the shock location. Although the flow vari- 
ables relax to the same values in the midst of the 
shock layer, the unresolved relaxation zone may still 
affect some important engineering variables, such as 
the radiative heating at the wall. The radiative emis- 
sion power behind the shock will depend strongly on 
the temperature and species densities, both varying 
rapidly in that region, and being very sensitive to 
the grid resolution. Therefore, a radiation code was 
used to compute the intensity along the line-of-sight 
in the stagnation region, and the heat flux at the 
wall: this computation was performed after the flow 
steady state was obtained, i.e. the flow and radiation 
were not coupled. After each adaption, the change 
in radiative heat flux was computed and compared: 
the results are shown in Table 2, for both the opti- 
cally thin case and the optically thick case. In the 
former case, the relative changes are quite large, and 
the values converge slowly. In the thick case, the 
absorption by the core of the shock layer tends to 
damp the perturbations: for that case, we see that 
the heat flux converges more rapidly towards a fi- 
nal value. Since the relative change is small (-1.8%) 
after the 3 rrf adaption, we considered that the res- 
olution was now sufficient. The comparison in ra- 
diative spectrum at the wall between the unadapted 
and final solution is shown in Figure 10. Most of the 
changes occur in the UV region. 
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It is important to remark that the computed radi- 
ation did not include the continuous spectrum, and 
therefore the variations in radiative heating at the 
wall may be under-estimated. Additionally, the den- 
sity is sufficiently large in this example that equi- 
librium radiation can be assumed: this considerably 
reduces the uncertainty in the computation of the 
radiative flux. For lower density and higher velocity 
cases, one must include thermal non-equilibrium ef- 
fects. The relaxation zone becomes then even more 
important to resolve accurately. 

It appears from this example that in practice it 
will be very difficult to compute a flow using an im- 
plicit method with a large CFL number, and that 
hydrodynamics-chemistry coupling effects will some- 
times limit the time step to CFL values below 1. Ad- 
ditionally, we may be required in practice to trans- 
form the grid, according to the solution obtained, in 
order to reach the desired accuracy: these adaptions 
need to be performed several times. It would seem 
therefore that a more efficient approach would com- 
bine the flow computation with the grid refinement. 
Indeed, there is a technique that can potentially lead 
to more efficient computations: using unstructured 
grids, the computational cells can be subdivided at 
will to give better accuracy in the regions that need 
it. Similarly, the subdivided cells can be regrouped 
in regions of low gradients, in order to keep the total 
number of cells within reasonable limits. Such a tech- 
nique would use a small number of cells to start with, 
and progressively refine them: most of the transients 
would then be computed using a small number of 
cells, leading to a more efficient procedure. 

VI: Conclusions and Recommendations 

We have not used here all the possible variations 
on the implicit schemes, neither have we exhausted 
the methods of coupling the chemistry and the fluid 
dynamics. We have however used techniques which 
are commonly used, and, we hope, demonstrated the 
trends for practical problems. We can draw several 
conclusions from this study: 

1. It is clear that on many problems of interest, 

the calculations cannot proceed with very large 

CFL numbers during the approach to steady 

state. Inevitably, for large numbers of species, 

the Block Tridiagonal methods cease to be ef- 

ficient in that regime. Only Scalar Tridiagonal 

methods, or even explicit methods remain effi- 
cient. 


2. It is clear that the Operator-Splitting approach, 
including the Explicit-Coupling between chem- 
istry and hydrodynamics, is at least as accurate 
as the Fully Coupled, and apparently more sta- 
ble for very stiff problems. Sub-iterations of the 
chemistry can further improve the accuracy and 
stability in the most severe cases of stiffness. 

3. The N s -Split method, at least in its present for- 
mulation, is too inaccurate for large time steps 
or strong concentration gradients. This disa- 
pointing result is not completely understood at 
the moment. It does not affect our conclusions, 
since the method is less performant than the 
Scalar method. This results should however be 
investigated further, since it may have implica- 
tions on other systems, such as a two- or three- 
fluid plasma, where the implicit treatment of the 
electron component gas dynamics is mandatory. 

4. Calculations of shock layers on fixed grids may 
not be sufficiently accurate if radiative phenom- 
ena or thermal non-equilibrium effects must be 
considered. In the example shown, several iter- 
ations at grid adaption were necessary. Other 
calculations on similar problems were also per- 
formed, that supported this conclusion. It ap- 
pears then that dynamical grid adapting should 
be performed during the course of the calcula- 
tion, for higher efficiency. 

Dynamic grid refinement could lead to even 
higher efficiencies if both the distribution and the 
overall number of grid points are allowed to vary. 
This can be done on structured as well as unstruc- 
tured grids. The construction of implicit schemes on 
unstructured grids would be quite complex. How- 
ever, we have concluded that this may not be a re- 
striction for many cases. An explicit algorithm will 
therefore be sufficient, and the technique is reduced 
to a sophisticated book-keeping problem. In addi- 
tion, the use of explicit, Operator-Split techniques 
allows us to take advantage of massively parallel (or 
mixed) computer architectures. This method will be 
investigated in the future. 

We have not mentioned another technique appli- 
cable for Operator-Split methods, when the chem- 
istry is very stiff. The chemistry (or other inter- 
nal process) can be rescaled, or ‘slowed-down 1 * * * * * 7 ar- 
tificially: this may have the effect of increasing the 
relaxation distances. However, we have made pre- 
liminary calculations that seem to indicate that in 
the very severe cases of stiffness, the changes are not 
perceptible. In addition, this procedure can be used 
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during the elimination of the transients, then the 
rescaling is progressively eliminated until a steady 
state with the proper time scale is obtained. If the 
rescaling is not eliminated, a false steady solution 
is obtained. The influence of the numerical proce- 
dure on the steady solution is also a serious ques- 
tion, discussed recently by Lafon k Yee [7], They 
show that for flows coupled with non-linear source 
terms, the steady state reached may depend on the 
path used to reach it. It is clear therefore that the er- 
rors induced by the numerical procedures can never 
be under-estimated, and that all users of CFD should 
proceed with extreme caution. 
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designation 

Treatment of hydrodynamics 

Chemistry/Convection coupling method 

FC/Block 
EC/Block 
OS/Block 
EC/Scalar 
EC/N, -Split 
EC-SI/Block 
EC-SIC/Block 

N x A r Block Tridiagonal (section II-A) 
N x N Block Tridiagonal (section II-A) 
N x N Block Tridiagonal (section II-A) 
Scalar Tridiagonal (section II-B) 
N s -split method (section II-C) 

N x N Block Tridiagonal 
N x N Block Tridiagonal 

Implicit (or Fully-) coupled 
Explicit (or Loosely-) coupled 
Operator-Split (or Fractional Step) 

Explicit coupling 
Explicit coupling 

Explicit with Sub-Iterations (or Sub-Cycling) 
Explicit, Sub-Iterations and Sub-Coupling 


Table 1: Designation of numerical methods and coupling methods used in this study. 


Grid Cycle 

Relative Change 

Relative Change 


in Surface Flux 

in Surface Flux 


Optically Thin Gas 

Optically Thick Gas 


[2000-8000 A] 

[1740-1750 A] 

i 

Non-adapted Grid - Adaption 1 

-25.1 % 

-11.0 % 

Adaption 1 - Adaption 2 

+27.1 % 

+5.2 % 

Adaption 2 - Adaption 3 

-3.5 % 

-1.8 % 


Table 2: Axisymmetric blunt body calculations - results of grid adaption study. 
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Required Memory (MW) Normalized CPU Time (CPUjmpiicit/CPUgxpiicit) 



Number of Species (N s ) 


Figure 1-a. Normalized CPU requirement for Block, Scalar and N s -Split methods 
(nonreacting). CPU is obtained as time per iteration per grid point, and normalized 
to CPU requirement for explicit method. 



Figure 1-b. Memory requirement versus number of species for same methods 
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Temperature (xlO 3 K) Temperature (xlO K) 



Mole fraction Mole fraction 


Temperature (xlO 3 K) Temperature (xlO 3 K) 


FC/Block 
EC/N S -Split 
EC/Scalar 


Distance from wall (cm) 

Figure 3-a. 1-D Shock (Case 2, P oe =10' 4 atm.). Results for CFL=2. 
Results of EC/Block and OS/Block methods are indistinguishable 
from FC/Block results and hence, are not shown. 


FC/Block 
OS/Block 
EC/N S -Split 


Distance from wall (cm) 

Figure 3-b. 1-D Shock (Case 2, P„=10‘ 4 atm.). Results for CFL=4. 
The EC/Scalar method failed for this case. 


Mole fraction Mole fraction 








Temperature (xlO 3 K) Temperature (xlO 3 K) 



Figure 4-a. 1-D Shock (Case 3, P_=10' 2 atm.). Results for CFL=2. 
The FC/Block method failed and the EC/N s -Split method gave 
nonphysical, oscillatory solution. 



Distance from wall (cm) 

Figure 4-b. Details of Case 3. FC/Block solution is obtained for a higher density 
(x 10) grid. Other solutions shown are for the standard grid but with subcycled 
chemistry or smaller time steps (CFL=0.2) 
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Mole fraction 





AE (J/m 3 ) AE (J/m 3 ) 



Figure 5-a. Quasi-ID Nozzle - energy residual history for the FC/Block method. 



Figure 5-b. Quasi-ID Nozzle - energy residual history for the EC/Block method. 
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AE (J/m 3 ) AE (J/m 3 ) 



Figure 5-c. Quasi-ID Nozzle - energy residual history for EC/N s -Split method. 



Figure 5-d. Quasi-ID Nozzle - energy residual history for EC/Scalar method. 
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Temperature (K) 



Distance along nozzle axis (m) 

Figure 6. Quasi-ID Nozzle - steady state axial distribution of atomic oxygen. 



0.0 0.2 0.4 0.6 0.8 1.0 


Time (xlO” 6 s) 

Figure 7. Quasi-ID Nozzle - Test of subiterated chemistry (EC/Block method). 
Initial state far from equilibrium and global time step (At) is 50 ns. Subcycled 
case (20 cycles) has an effective time step (SO of 25 ns. The temperature is 
recomputed only at the end of the integration over the global time step 50 ns. 
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Molar Concentration (moles/m ) 





"3 

Distance from body surface (xlO" m) 


Figure 8-a. Axisymmetric Blunt Body - Temperature 
profile along the stagnation line (non-adapted grid). 



Distance from body surface (xlO 3 m) 

Figure 8-b. Axisymmetric Blunt Body - Electron mole fraction 
profile along the stagnation line (non-adapted grid). 


18 



Mole Fraction of e’ (xlO 3 ) Temperature (xlO 3 K) 



Distance from body surface (xlO m) 


Figure 9-a. Axisymmetric Blunt Body - Temperature profile along 
the stagnation line for non-adapted and adapted (final pass) grid. 
Note the significant change in shock location and peak temperature 



Figure 9-b. Axisymmetric Blunt Body - Electron mole fraction profile 
along stagnation line for non-adapted and adapted (final pass) grids. 
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Wavelength (A) 

Figure 10. Axisymmetric Blunt Body - Comparison of radiation spectra 
incident on the wall (stagnation line of sight). The change is = 35% in 
the UV(2500 A) region. 
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Abstract: 

This paper describes the computational work be- 
ing performed at Ames on the simulation of the 16” 
Shock Tunnel facility. The paper describes the ap- 
proach used and shows some preliminary results for 
various flow transients. In particular, we describe the 
numerical problems encountered during the compu- 
tation of these flows, and the methods used to resolve 
them. We also discuss the validity of some approx- 
imations used, notably concerning the reduction of 
the problem into problems of smaller dimensionality, 
or smaller size. We show how quasi- ID simulations 
can be used to help design experiments, or to better 
understanding the characteristics of the facility. An 
application to the design of a non-intrusive diagnos- 
tic is shown. The multi-dimensional flow transients 
computed include the shock reflection at the end of 
the driven tube, the shock propagation down the noz- 
zle, and the breaking of the main diaphragm. The 
interaction between separate flow events will also be 
discussed. 


I. Introduction 

The Ames 16” facility, shown schematically in 
Figure 1, can be considered as a typical example of a 
shock-tunnel for hypersonic flows. The shock tunnel 
is about 70 meters long, composed of a driver sec- 
tion (17” diameter) filled with a light combustible 
mixture, and a long driven tube (filled with the test 
gas at low pressure). At the end of the tunnel is a 
supersonic nozzle 6 meters long (1 m diameter at the 
exit pl£ne) and finally a test section. The complete 
operation of the facility typically results in test times 
of the order of 5 to 20 milliseconds: the time for the 
main shock to propagate from the main diaphragm 
to the end of the driven tube is of the order of 7 msec. 
After partial reflection at the end of the driven tube, 
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the shock interacts with the incoming contact discon- 
tinuity (CD) which separates the test gas from the 
driver gas. For ’tailored’ conditions, there is no wave 
reflected back from this interaction, and the flow be- 
tween the CD and the nozzle is uniform and steady. 
More detailed descriptions of a shock tunnel opera- 
tion can be found for example in ref. [1]. The test 
time is measured from the time the flow conditions 
in the stagnation region become steady, until the CD 
reaches the end of the driven tube, and the nozzle 
flow becomes contaminated with the driver gas. A 
flow ’quality 5 will be represented by the steadiness of 
the stagnated flow, the low level of contamination of 
the test gas by driver gas or impurities, as well as 
the peak conditions (pressure, enthalpy) attainable 
by the facility. 

Numerical simulations are required to better un- 
derstand the test flow conditions (temperature, mass 
fraction profiles, time dependence, etc..) and sup- 
plement the experimental knowledge: they are used 
also to improve the design of the experiments or to 
improve the tunnel operation characteristics. In or- 
der to satisfy these objectives, several aspects of the 
tunnel operation need to be simulated. Modeling 
the entire facility from the driver to the test section 
requires us to grid a physical length of about 50 me- 
ters; yet, important flow features such as shock and 
contact discontinuities (CD) should be resolved^vith 
good accuracy, i.e. a few mm in the axial dire<?flE7n. 
The spatial stiffness is then of the order of 1 0 4 ; this 
is a conservative estimate, since the boundary layers 
also need to be resolved in some cases, requiring grid 
spacing as low as tens of micrometers in the radial 
direction. The modeling of the shock-tunnel facil- 
ity falls into the general category of multiple scale 
problems, and calls for appropriate special numeri- 
cal methods which we will hint at in this paper. 

The complete operation of the hypersonic facility 
also involves a large number of different physical pro- 
cesses at work, some of which are not well understood 
or are very difficult to compute. The combustion 
process in the driver gas requires a good description 




of the energy deposition and flame propagation; the 
main diaphragm rupture would require us to under- 
stand the material deformation up to the plasticity 
limit. The penetration of the jet of driver gas into 
the driven tube is a problem of 3D, turbulent, multi- 
scale mixing. The temperatures are sufficiently high 
that some wall ablation takes place and contaminates 
the flow. Finally, radiative effects may need to be 
considered, and low density, thermal nonequilibrium 
gas models may be required in the nozzle expansion. 
A multiplicity of CFD tools must therefore be made 
available to study the effects of each of these phe- 
nomena. 

The problem is compounded by the lack of cru- 
cial experimental data: such phenomena such as di- 
aphragm rupture are difficult to observe, and only 
rough estimates of the process time and energy scales 
can be made available and approximately correlated 
with tjie experimental data. Ablation of the tun- 
nel material is also a very complex physical process, 
which may depend on the microstructure of the ma- 
terial itself. 

Solving this problem will require us to develop 
and test appropriate numerical techniques, decom- 
pose and reduce the problem accordingly, and de- 
velop and test various physical models. Simplifica- 
tion of the problem can be accomplished first by ’di- 
mensional reduction 1 , i.e. solving the problem in in 
a reduced number of dimensions. Quasi-ID simula- 
tions have low computational time and memory re- 
quirements, and therefore can be used for a large 
number of computations, such as those required for 
design and sensitivity studies. It will be important 
to verify the accuracy of these simulations, and to 
estimate the quality of information which can be ob- 
tained from them. 

In a typical cartesian tradition, decomposition 
of the problem into several independent sections is 
usually attempted. This we regard as ’causality re- 
duction 1 . This approach requires several approxima- 
tions and simplifications, and its validity needs to be 
more clearly justified. For example, the time evolu- 
tion of the flow in the test section clearly requires 
us to accurately compute the unsteady flow in the 
stagnation region, including gas/wall interactions, 
shock/boundary layer interactions. The shock/CD 
interaction is also very important, and this requires 
us to know very well the shape of the CD. Going 
back further, the shape of the CD will be influenced 
by the earlier process of flow establishment from the 
main diaphragm rupture and thereon. But the de- 


tails of the diaphragm rupture will depend on the 
pressure rise in the driver gas, and this requires us 
to model first the turbulent flame propagation in the 
driver tube. The logical sequence of events, i.e. from 
driver tube combustion to flow expansion into the 
test chamber, is unfortunately the most difficult to 
compute, since it requires us to be able to model the 
most complex phenomena at first. Therefore we have 
to sacrifice some degree of detail in order to obtain 
practical answers in reasonable time. By causally re- 
ducing the problem, we loose accuracy but gain in 
efficiency: the method can still be used to examine 
the important physical effects (i.e. sensitivity stud- 
ies) in an efficient way. 

In this paper, we will discuss these methods and 
the preliminary results. The following questions 
must be addressed: 

• what improvements in numerical techniques are 
required. 

• how accurate are the results from reduced (ID) 
models. 

• how significant are the causal influences of var- 
ious sections of the tunnel. 

• what physical phenomena are important and 
how to model them. 

The last item on our list will not be discussed 
here. In this paper, we will concern ourselves primar- 
ily with the numerical techniques, and the evalua- 
tion of the approximations commonly made in shock- 
tunnel simulation. The obvious intent of the devel- 
opment of this numerical capability is to provide an 
array of numerical tools for better understanding of 
the facility operation, and the design of new test con- 
ditions and diagnostic procedures. The combined fo- 
cus of both experimental and numerical tools leads to 
superior measurement capability, and is an approach 
practiced at NASA- Ames [2], 

II. Numerical Techniques 

We have already mentioned that the spatial and 
time stiffness can be quite large; fortunately the CFD 
community has developed an arsenal of techniques 
which can potentially be used to reduce the severity 
of the problem. First, we observe that the flow dis- 
continuities are few and generally well localized: it 
is therefore unnecessary to simultaneously carry the 
same resolution requirements in all regions of the fa- 
cility. We can rely on several techniques of dynamical 
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grid refinement to carry this task. Second, the prob- 
lem can be approximately separated in three parts. 
The combustion process in the driver gas proceeds 
independently of the driven tube, and can be com- 
puted separately if necessary (throughout this paper, 
it will be simply modeled). Since the flow is mostly 
supersonic in the nozzle, there is no backward influ- 
ence of the flow from the nozzle to the tubes. It is 
therefore adequate to compute at first the unsteady 
flow in the driveT and driven sections, including the 
nozzle and up to the throat; the flow solution at the 
exit plane can be stored according to the required 
spatial and time accuracy, and used as unsteady in- 
flow condition for the remainder of the facility, i.e. 
nozzle and test section. This allows us also to switch 
between different physical models (e.g. nonequilib- 
rium thermodynamics) and different dimensionalities 
(3D versus 2D axi-symmetric or 2D versus ID). 

Another important question concerns the choice 
of numerical method: although we are interested 

in the unsteady flow motion, the flow time scales 
that need to be resolved may be large compared to 
the time scale from the stability limit of an explicit 
scheme. This is true especially for a clustered vis- 
cous grid. An implicit method, if time-accurate, can 
potentially be used for this problem as well. How- 
ever, there is a huge cost associated with the inver- 
sion of the block-tridiagonal matrices from the LHS 
of the system of equations, especially when multi- 
ple species are present. Our experience showed that 
even when only two species are considered (driver 
gas/driven gas), the implicit method needs to be ex- 
ecuted at a CFL number > 5 for better performance 
than the explicit method. Due to the severe nature of 
the problem (pressure ratio, shock speed), this CFL 
value is generally too large for these transient flows. 
However, the viscous terms alone can be treated im- 
plicitely, and show greater stability than the convec- 
tive terms: our method therefore performs an Opera- 
tor Splitting between convection, viscous dissipation 
and chemistry. The viscous terms are treated im- 
plicitely, and the inversion of the 3 x 3 (in 2D) block 
tridiagonal matrices is sufficiently fast to be justified. 
The chemistry is also treated point-implicitely, and 
the algorithm is also very fast. When the chemistry 
is extremely stiff however, linearization errors will 
reduce the stability of the point-implicit algorithm. 
For example, let us consider air (79% N 2 , 21% O 2 ) 
suddenly raised to 100 atm and 6000 °K: at con- 
stant volume, the system reaches equilibrium within 
a fraction of a microsecond. Attempting to solve 
implicitely the chemistry in one iteration and with a 
time step larger than ~ 2 10 _s sec would lead to very 


unstable and unphysical solutions. The initial evo- 
lution from the highly nonequilibrium state needs to 
be time-resolved more accurately: our algorithm sub- 
iterates (more precisely, ’sub-cycles’) the chemistry, 
i.e. integrates over a given time interval At with vari- 
able steps St . The sub-step is initially estimated from 
the rate of change of the species, then is stretched in 
a fixed proportion for the subsequent sub-iterations: 
5t n + l = (1 + a)6t n } with a > 0. Each sub-iteration 
is solved point-implicitely. The change in tempera- 
ture is also estimated at each sub-iteration from the 
change in chemical energy. Figure 2 shows a sample 
computation (a = 0.5) for this highly nonequilib- 
rium case. The ’exact ’ time evolution is also shown. 
The integration intervals are At = 5 10“ 7 seconds. 
One can see that the variable step solution is off in 
the equilibrium region until the end of the first inte- 
gration interval, at which time the thermodynamic 
state of the gas is re- analyzed and a more accurate 
calculation of the temperature is performed. The 
sub-iterated solution then quickly matches the exact 
solution. 

To maintain good accuracy, the relative changes 
in temperature estimated during the chemistry 
should be limited. A chemical time scale is there- 
fore defined as the time during which the chemical 
reactions induce a relative change in temperature of 
~ 5%t . The effect of chemistry on the flow dy- 
namics occurs principally through the temperature 
change, and this time step limit provides us with 
a criterion for accurate coupling of the chemistry 
and flow dynamics. If there is no monitoring of the 
chemical time scale and no enforcement of this time 
step limit, an instability generally develops rapidly, 
especially in cases of stiff 1 chemistry; this instabil- 
ity can occur for both a fully-coupled approach or 
the Operator-Splitting approach. The case of very 
stiff chemistry would seem therefore to be very in- 
efficient: this corresponds basically to a very poor 
resolution of the chemical relaxation distance, such 
as the one behind the primary shock. However, it 
is possible to somewhat rescale the chemical time 
scale without changing the solution too much. Let 
us suppose that we want to resolve the flow dynamics 
on the order of a convective time-scale At conin but 
the chemistry is so stiff that accuracy and stability 
considerations restrict us to (for example) a global 
time scale ~ 0.01A< conv . We can artificially rescale 
the chemistry by 10, and the stability limits lead us 

tThis number may be varied, according to the desired 
accuracy and/or the stability. It is our experience that in 
most cases, the maximum relative change allowed should be 
less than 10%. 
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to a time scale 0.1A* conu , which we can better af- 
ford. This rescaling can be done locally, depending 
on the stiffness ratio: this will affect slightly the so- 
lution, for example by over-estimating the chemical 
relaxation length behind a shock, but the difference 
may be practically insignificant for very stiff cases. 
The technique is easily implemented by restricting 
the integration interval (it becomes 5 10~ 8 sec in our 
example): the insert in Figure 1 shows the results 
of this rescaling approach. We see that the original 
transients are still reproduced, and the solution still 
converge to the exact solution. This sub-iterated al- 
gorithm, with or without the rescaling option, is used 
in all our calculations. 

The spatial stiffness can be solved by various tech- 
niques of dynamic grid adaption. In the driver tube, 
strong gradients and discontinuities exist only dur- 
ing the combustion process. After rupture of the 
main diaphragm, only weak gradients subside, and 
this section of the facility does not require high reso- 
lution. The spatial accuracy is required for the prop- 
agating shock and CD, down the driven tube. One 
possible technique is to dynamically compress and 
stretch the grid points to accumulate them in the 
required regions. This technique can work well in 
one dimension and is easy to implement: care must 
be taken however to prevent sudden jumps in spac- 
ings, or accuracy will be lost. The technique may be- 
come more problematic in two or more dimensions, 
as some flow features (including the CD) may have 
more complex shapes: the grid can become distorted, 
and can affect the solution. Another technique, pio- 
neered by a group at Livermore [3], consists of adding 
smaller scale grids in regions of interest. These grids 
can be exact sub-scale replicas of some regions of 
the coarse, background grid: the flow variables can 
be transferred between grids in an exactly conser- 
vative way. This way, the grids are never distorted 
and their motion can be computed very quickly, with 
minimal overhead. An example is shown in Figure 3: 
the propagation of the primary shock and CD down 
the driyen tube is computed in a one-dimensional 
model. A first calculation used a constant spacing 
grid, with about 1000 points for the whole tunnel. 
A second calculation used a background grid for the 
whole tunnel with 250 points, and a high resolution 
grid super-imposed on it: the latter moves along with 
the shock, with a sliding motion. The sub-scaled grid 
had a size of 800 grid points, and the scale ratio be- 
tween the two grids was 20 to 1. Flow values within 
the background cells are computed by volume aver- 
aging of the sub-scaled cells present, if any. While 
both calculations used approximately the same num- 


ber of points, it is clear that the second method has 
a much higher local resolution, and the sharpness of 
the CD is dramatically improved. 

Another method, which we will be testing in the 
future, uses non conservation of the number of grid 
points. Unstructured grids can be manipulated to 
accumulate points in the regions of interest, either 
by grid displacement or by dynamic subdivision. 

The singular axis was found to be another recur- 
rent problem during the computation of flow tran- 
sients: the simulation of the shock reflection at the 
end of the driven tube, for example, initially showed 
a strong conical shock structure near the axis with 
the apex of the cone leading the overall structure (see 
Figure 4). This peculiar formation can also be seen 
for example in the results of P. Jacobs [4]. This struc- 
ture is possible only if a very intense and high veloc- 
ity jet of gas is produced and maintained on the axis: 
this is a highly singular and unphysical behaviour. 
Close examination of the numerical results indeed 
showed an excessively large axial velocity component 
in a single cell close to the axis. Because this high 
velocity jet was present in one cell only, and did not 
show signs of diffusion, we were convinced that it 
is the result of a numerical error. A similar phe- 
nomenon can be observed also in the propagation of 
the main shock down the nozzle [5], and can be seen 
at the axis or in some cases near the walls. This phe- 
nomenon was observed only for axi-symmetric flows, 
and when a second-order accurate scheme (with min- 
imal dispersion) was used. It was finally related to 
the aspect ratio of the grid cells: indeed, the problem 
disappeared when the aspect ratio of the grid cells 
was adjusted to lower the radial gradients. If the 
spacing is such that Ax << Ar, the flow features 
are smooth and accurate. If Ax > Ar an instability 
may develop in some regions of the flow. We assume 
the problem comes from the axi-symmetric pressure 
correction term, which is not part of the monotonic 
(TVD) fluxes and acts as a non-conservative momen- 
tum source term. This pressure correction term is a 
result of the formulation of the Euler equations for 
an axi-symmetric problem, and can be easily visu- 
alized with the following finite- volume description. 
Let us consider a cell in a cylindrical geometry, de- 
scribed schematically in Figure 5. As one approaches 
the axis, the ratio of cell side surface to cell volume 
behaves as 1/Ax for the axial direction, while it is 
exactly 0 for the radial direction. There is a contribu- 
tion to the radial momentum density from the pres- 
sure on the sides of the cells in the azimuthal direc- 
tion. This contribution is proportional to 1/Ar. It is 
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clear that when the cells are clustered near the axis, 
the contribution from the non-conservative momen- 
tum source term may become dominant, and there- 
fore there is no guarantee that a non-oscillatory solu- 
tion can still be obtained. The problem usually dis- 
appears when the grid is relatively coarse in the ra- 
dial direction, near the axis, in which case the source 
term ceases to be the dominant contribution. This 
observation also applies to other cases, such as the 
computation of shock layers on blunt axi-symmetric 
bodies. In most unsteady cases or for steady flows, 
the radial gradients will be weak near the axis. One 
particular exception would be the propagation of a 
shock ring towards the axis: this is such the case 
in the reflection of the primary shock in the tunnel, 
from the upper section of the end wall. Figure 6 
shows a time sequence for the same case as Figure 4, 
but for a grid coarsened near the axis: we see that 
the instability is effectively removed. A similar im- 
provement can be obtained by reducing the spatial 
order of accuracy in that region. 

The general effect of the grid structure or aspect 
ratio on the flow solution should never be under- 
estimated, especially in the severe cases we are con- 
cerned with. Although we may see unstructured 
grids as a potentially useful tool for the flows of in- 
terest here, their effect on the flow solution will have 
to be carefully evaluated. 

The flow expansion in the supersonic nozzle is one 
of the easier tasks to perform: the transient flow sim- 
ulation is required to examine the steady flow estab- 
lishment time in the test section, and to verify that 
the flow does not choke during that time. This is 
particularly relevant for large objects or for massive 
gas (fuel) injection. Since the nozzle is evacuated to 
very low pressure before the rupture of the secondary 
diaphragm, the flow propagating down the nozzle is 
preceded by a high velocity jet. The mean free path 
in that region can be quite large (~ 3 mm for P = 100 
milliTorr), and strictly speaking the Navier-Stokes 
equations cease to be valid in this low-density re- 
gion prior to the shock. By attempting nevertheless 
to solve the flow dynamics with the Navier-Stokes 
equations, we are experiencing strong viscous effects 
which operate on a very short (~ 10 pico-second) 
time scale. Although the viscous fluxes are computed 
implicitely, it is necessary for stability reasons to ar- 
tificially reduce the strength of these viscous effects 
in that region. This is easily accomplished by us- 
ing a numerical switch that effectively and smoothly 
removes the gradients on a scale smaller or compa- 
rable to the mean-free path during the calculation 


of the viscous fluxes. This switch will effectively re- 
duce the spreading of the shock into the low density 
gas. The shock thickness will therefore be under- 
estimated. The cutoff length scale can be adjusted 
arbitrarily: in Figure 7 we show the effect of the cut- 
off choice (0.1 versus 10 mean free paths) on the solu- 
tion, for two cases of nozzle pressures. The change is 
dramatic for the first case (100 /zTorr) and unnotice- 
able for the second (100 milliTorr). Nevertheless, for 
the latter case, which is also the case we will be using 
throughout the remainder, the viscous time scale is 
reduced by a factor of 5 by choosing the higher (10 
m.f.p.) cutoff value. 

Using all the modified techniques now at our dis- 
position, we are able to sucessfully compute the tran- 
sient flows in the nozzle and driven tunnel at the 
rupture of the secondary diaphragm (see Figure 6). 
Different grids were used for the driven tube and 
nozzle regions, and no subscaled grids were used for 
this particular case: both regions were coupled at 
each iteration. The results of the computation for 
the nozzle flow (Figure 8) are in good qualitative 
agreement with an experimental schlieren of a simi- 
lar problem shown in Figure 9, taken from ref. [6]. 
Both show the curved leading shock, and a complex 
structure behind it, dominated by a Mach disk, it- 
self supported by oblique shocks emanating from the 
nozzle walls. Although the initial conditions for these 
two flows are very different, the similarity between 
the two structures lead us to conclude that they are 
examples of a general pattern. Details of the con- 
ditions (pressures, geometry, etc..) may change the 
relative strength and position of the shocks, with- 
out modifying the overall configuration. A snaphost 
taken at later times would show that near the exit 
plane of the nozzle, the primary shock straightens 
out and becomes normal. If one was to perform an 
unsteady computation in the test chamber assuming 
a normal shock at the inflow, the calculation would 
be in error by leaving out the complex shock struc- 
ture which propagates immediately behind the lead- 
ing shock. This is shown for example in Figure 10, 
where a cone has been used as testing body. 

Ill, Dimensional Reduction 

When test times are large compared to the tran- 
sients, it seems appropriate to compute the nozzle 
flow for the steady state, without having to perform 
the calculation time-accurately for tens of millisec- 
onds. This allows us to use many numerical tech- 
niques to make this computation more efficient. We 
have done several such computations, solving for the 
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full Navier-Stokes equations, but for limited regions 
at a time. The computation proceeds on a subset 
of the whole nozzle grid, until convergence is ob- 
tained. The subset then moves further down the 
nozzle; the procedure is similar to the displacement 
of a ’window’ along the flow direction. This proce- 
dure is half-way between a global calculation and the 
PNS method; although less efficient than the PNS 
method, it is more flexible and allows us to correctly 
treat embedded recirculation zones or other unex- 
pected subsonic regions which may occur. It was 
observed that the convergence rate and the compu- 
tation time are dominated by the chemistry in the 
boundary layer, and the downstream propagation of 
chemical changes in the low velocity sublayer. Al- 
though the calculation does not display severe in- 
stabilities such as may be the case for the transient 
flows, it is still desirable to coarsen the grid in the 
radial direction near the axis, to avoid spurious fluc- 
tuations in pressure. 

The steady nozzle flow computations are usually 
done to obtain estimates of the flow conditions at the 
exit plane, in particular temperature, Mach number, 
pressure and species concentrations. Since it is gen- 
erally assumed that the flow is nearly uniform at the 
exit plane, it seems that simple one-dimensional com- 
putations should be adequate. We have performed 
such computations using a quasi- ID version of our 
code, and compared the results. The effect of the 
boundary layer is expected to be the major contri- 
bution to any potential disagreements between the 
quasi- ID solution and the 2D axi-symmetric solu- 
tion. This effect is taken into account by effectively 
modifying the nozzle profile and its area. By con- 
sidering the inviscid core only, the quasi- ID method 
can approximately take into account the constricting 
effect of the boundary layer by assuming a new noz- 
zle shape for which the effective radius is a constant 
fraction of the real nozzle radius. The dependence 
of some flow quantities on the effective nozzle area is 
shown in Figure 11. It is clearly apparent that the 
static pressure has the highest sensitivity to the ef- 
fective nozzle area, and therefore to the boundary 
layer thickness. Temperature and species concen- 
trations are less sensitive to the area variation, and 
additional uncertainties about the chemical rates or 
contaminants are likely to cloud the issues: the static 
pressure is therefore an important variable to mea- 
sure at the nozzle exit. The comparison with the 2D 
axi-symmetric results is shown in Figure 12, where 
again the sensitivity of the static pressure is clearly 
demonstrated. The best agreement is obtained for 
an effective nozzle radius of 87% the real radius (i.e. 


a boundary layer thickness approximately 13% the 
nozzle radius). Notice also that the static pressure 
is the only variable that still displays significant ra- 
dial fluctuations within the inviscid core, at the exit 
plane. A multi-point measurement of this variable is 
therefore doubly informative, as far as code valida- 
tion is concerned. In Figure 13 we show the radial 
profile of Temperature and NO concentration at sta- 
tion code-named N3 (2.37 meters downstream of the 
throat) and the exit plane. All the species profiles 
are relatively unchanged between N3 and the nozzle 
exit, except atomic nitrogen N t which is present in 
very small amounts. Owing to the high reactivity of 
N , this result is not surprising: most of the chem- 
istry is frozen before station N3, and this is especially 
true of NO. The temperature variation between N3 
and the exit plane is therefore due to hydrodynami- 
cal effects only. 

NO is a relatively important specie to measure 
(it effectively ties up a significant amount of oxygen 
and has a noticeable effect on ignition delays), and it 
can be done easily with a laser absorption method. A 
computational study of this diagnostic technique can 
also be done to help design the experiment. Th^ com- 
puted high resolution spectrum of the NO(j) band 
system is shown in Figure 14, with and without the 
boundary layer. The intensity at peak absorption of 
the (0,0) band changes by a factor of ~ 3 when the 
boundary layer is removed. A precise measurement 
of the core flow must avoid the uncertainty caused by 
the boundary layer: this can be simply done by pro- 
viding an optical wave guide that protects the beam 
from the boundary layer. Although this protrusion 
is likely to disrupt the flow, the effect may be small 
enough or irrelevant at the exit plane. We can now 
look at the effect of core flow variation by comparing 
the intensity for the 2D solution (without BL) and 
a quasi- ID solution. The result is plotted in Fig- 
ure 15-a as the ratio of intensities between the two 
cases. Surprisingly, the error is less than 10% for the 
most sensitive lines. At the exit plane (Figure 15-b), 
the results are even better, with deviations less than 
2%! This seems to indicate that this measurement 
can yield useful information on the core flow, with- 
out having to worry too much about details of the 
radial fluctuations within the inviscid core. 

TV: Causality Reduction 

In the previous section we have seen that dimen- 
sional reduction of the problem can still yield some 
useful results, and quasi-ID calculations can help de- 
sign and understand the facility operation and exper- 
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imental diagnostic techniques. Some flow features 
are lost in this reduction process, for example the 
shock structure that propagates down the nozzle. If 
one were to compute the transient flow inside the test 
section, one would require a prior computation of the 
transient nozzle flow. This is an example of causal 
interaction between two different regions of the fa- 
cility. The computation of the transient nozzle flow 
can be done by assuming constant stagnation con- 
ditions, or starting from a uniform shock impinge- 
ment on the end wall of the driven tube. For later 
times, the propagation of the reflected shock in the 
driven tube must be followed accurately: this entails 
the interaction of the shock with the bounday layer 
and with the contact discontinuity. If this numeri- 
cal simulation is done with sufficient detail, we can 
gather information on the effect of boundary layer 
and CD shape on pressure wave generation and flow 
contamination. The shape of the contact discontinu- 
ity may however depend on its early history, i.e. on 
the main diaphragm rupture itself. This is another 
example of causal interaction, which forces to model 
the diaphragm rupture and subsequent evolution, in 
order to at least estimate the effect on CD shape. 
This is another difficult problem which will require 
a lot of effort: the calculation done so far is only a 
first attempt at solving the problem, which helped 
us identify the areas where further improvement is 
necessary. 

Although the opening and petaling of the steel di- 
aphragm is a fundamentally 3-dimensional process, 
it is not required at this stage to include this compli- 
cation. The problem is therefore reduced to a 2D 
axi-symmetric problem, and the diaphragm open- 
ing is reduced to a case of time varying boundary 
condition. The driver and driven tubes are grid- 
ded as two separate regions: between them, the 

boundary condition is set as a reflecting wall for 
some of the grid cells, and a patching condition for 
the grid cells within the opening. The distribution 
of patched/reflecting boundary points changes with 
time; in this first attempt, the grids are assumed 
fixed in time, i.e. the opening proceeds in steps, one 
grid point after another. If the grid had sufficent res- 
olution in the radial direction, this would be a good 
approximation to a continuous process. Some effects 
are ignored in this approximation. First, the physical 
boundary (steel) between the two regions is changing 
in time, due to the distortion of the diaphragm (the 
petaling); its dynamics should be modeled as well, 
and this will require a major effort in the future. We 
also assumed that the diaphragm is initially vertical, 
while in reality it is to some extent assuming a hemi- 


spherical shape, due to the pressure rise in the driver 
tube. Since we are only interested at this point in 
testing the numerical capabilities, this is not a con- 
cern. 

In Figure 16 we show the temperature contours 
of the flow transient, taken at three different times 
after the start of the opening of the main diaphragm. 
At 30 //seconds, the primary shock is just past the 
contracting section of the driven tube, immediately 
followed by the contact discontinuity and a hot layer 
of driver gas. The opening at that time is still small 
(about l/Q ih of the driven tube diameter). Notice 
that between the hot, sheared gas layer and the open- 
ing is a weak shock (blue contours in the figure) is 
emanating from the diaphragm opening. This is the 
result of the step in diaphragm opening: as one more 
cell of the driver section is put in contact with the 
driven tube, the sudden change in local boundary 
condition creates this weak shock. This purely nu- 
merical effect was difficult to estimate a priori , and 
is an interesting observation in itself: it forces us to 
reconsider the technique for future computations, in 
order to have a smoother, continuous opening. At a 
later time (54 // seconds), the primary shock has trav- 
elled further down the driven tube and significantly 
straightened. Oblique shocks reflected from the walls 
of the constricting section interact near the axis to 
form a strong mach disk. The shear layer has at- 
tached itself to the walls, and the driver gas behind 
the primary shock has developed a concave shape 
near the axis. At 70 //seconds, the primary shock 
is completely planar, the mach disk has shrunk, and 
a very complex flow structure follows the primary 
shock and CD. Since the stepwise diaphragm open- 
ing has produced spurious shocks, it is difficult to 
identify the real physical effects. Presumably, this 
structure is actually simpler (at least in this 2D ap- 
proximation). We can safely assume that a few main 
conclusions can be drawn from this preliminary re- 
sult: 

• the primary shock becomes planar very rapidly 
(~ 60/xseconds). 

• a complex and unsteady flow structure, domi- 
nated by a mach disk, is formed behind the CD 

• the CD itself develops a complex shape. 

The last item can be more clearly demonstrated in 
Figure 17, which shows the CD and the main shock, 
for the same time sequence. It is clear from Figure 
17-c that the CD is definitely non-planar; indeed, it is 
likely that from that moment on, the CD evolution 
will be dominated by Rayleigh-Taylor instabilities, 
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and will not recover an ideal, planar shape. 

The remainder of the simulation, i.e. the full 
propagation of this structure down the whole length 
of the driven tube, would take a considerable amount 
of CPU, and has not yet been attempted. Future 
work will focus again on the diaphragm opening and 
on a numerical technique that more accurately de- 
scribe it. It is important also to point out that this 
calculation is difficult for an additional reason: as 
the diaphragm opens, there is a strong jet of gas 
that flows through a small opening in a ideally re- 
flecting surface. This jet has a tendency to entrain 
some flow around the opening, producing an artifi- 
cial ’cavitation’ near the surface. The same numer- 
ical phenomenon can be also observed occasionally 
near the base wall of a wedge in a hypersonic flow, 
for example. This numerical instability can be re- 
moved by enhancing the numerical diffusion in that 
region- This was done so in our case. This will also 
affect the diffusion of the driver gas into the driven 
tube, and therefore the extent of mixing and shape of 
the contact discontinuity. This transient flow must 
therefore be recomputed with special care. 

Another example of causality relation concerns 
the influence of the combustion process in the driver 
and the flow in the driven tube. Pressure fluctua- 
tions are experimentally observed in the driver tube, 
which seems to indicate that the combustion is not 
a uniform process, and that some pressure waves are 
bouncing back and forth in the driver tube. These 
waves can travel down the driven tube after opening, 
and influence the pressure field in the stagnation re- 
gion. This effect is demonstrated in Figure 18, where 
the pressure histories at two locations are compared 
with the computations: Figure 18-a shows the pres- 
sure history a few centimeters upstream of the main 
diaphragm, in the driver tube. Figure 18-b shows 
similar profiles near the end wall of the driven tube. 
Pressure fluctuations with a sine shape have been 
superimposed in the driver tube at the moment of 
diaphragm rupture; the amplitude of this fluctuation 
was chosen to match the experimental observation in 
the driver tube. These quasi- ID computations show 
a very good agreement with the experimental traces 
when the fluctuations (of the right phase) are super- 
imposed. The agreement is notably much better than 
without these fluctuations (Figure 18-b). It is clear 
that the combustion process in the driver tube has a 
strong influence on the flow conditions in the stagna- 
tion region, and should also be modeled and better 
understood. These calculations were performed by 
G. Wilson at NAS A- Ames, using a quasi- ID code: 


additional details on the modeling will be presented 
in the future [7]. 

These results are another example of the useful- 
ness of dimensional reduction: it also shows that with 
some clever modeling, one can causally reduce the 
problem and retain an accurate description of the 
system. This should also help in understanding the 
conditions for better reproduction of test runs, and 
for more uniform flow conditions. This is another 
difficult task, which involves the modeling of energy 
deposition, initiation of deflagration and flame prop- 
agation. An effort in that direction is planned. 

V; Conclusions and Future Plans 

From these preliminary results, we can draw the 
following conclusions: 

• The flow conditions are very severe and put 
enormous strain on the accuracy and stability 
of the current numerical techniques. Further re- 
search into the improvement of the numerical 
techniques is desirable. 

• The validation of the CFD capabilities will re- 
quire some difficult measurements, including 
high resolution video recording of the tran- 
sient processes, such as secondary and main di- 
aphragm rupture. 

• The numerical modeling of the facility, even 
with simplistic assumptions (quasi- ID), can 
greatly benefit the design of experiments, diag- 
nostic procedures, new test conditions, and un- 
derstanding of the tunnel performance. 

• The computation of separate regions of the facil- 
ity treated as independent is an approximation 
at best. The influence of remote and past events 
on the overall flow structures is not negligible, 
and must be estimated. The separation of the 
facility flow into separate regions can be used for 
the estimation of various physical phenomena. 

The most immediate challenge still concerns the 
opening of the main diaphragm. The preliminary re- 
sults shown here suffer from an inadequate treatment 
of the unsteady boundary condition, and must be 
recomputed more accurately. Other challenges that 
await us concern the importance of various physical 
effects, such as shock-BL interactions at the end of 
the driven tube, and wall chemistry and ablation and 
resulting flow contamination. These three topics will 
be addressed in the future. 
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Figure 1: Schematic (not to scale) of Ames 16” Shock Tunnel facility. All dimensions are in cm. 



Figure 2: Computation of stiff chemistry. The open symbols are obtained when the implicit chemistry 
solver is sub-iterated with an increasingly large sub-step. The temperature change is also estimated at each 
sub-step. The integration proceeds until At = 5 10" 7 seconds, which is an assumed global time step. Deviations 
from the ’exact’ solution at times greater than 10" 7 seconds are due to an error in the estimate of temperature. 
As the first integration is terminated, a more accurate estimate of the temperature (including variations of the 
specific heat) is made, and the solution quickly converges to the correct one. The insert shows the effect of 
rescaling of the chemical time scale: here the integration proceeds only to 10% of the global time step. The 
integration during the second global step also quickly converges to the exact solution. 
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Figure 3: Grid subscaling performance. A constant fixed grid (1000 points) is compared with the results 
from a high density grid (800 points) moving along with the shock. This final profile, at the end of teh driven 
tube, shows the increase in resolution of the contact discontinuity. The high density (’subscaled’) grid is a subset 
of the ’background’ grid (250 points), and moves by steps equivalent to one cell spacing of the background grid. 
The transfer of information between the two grids proceeds by volume averaging, and is fully conservative. 





Figure 4: Driven/Nozzle flow transient at 80 
fisec after shock arrival. The shock comes from 
the left, instantaneously ruptures the secondary di- 
aphragm, propagates down the nozzle and is partially 
reflected back into the driven tube. Cells are clus- 
tered near the axis, the calculation i<jf for inviscid, 
non-reacting flow, spatially 2 n<i -order. The conical 
structure is believed to be an artifact from the nu- 
merical method. 
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Figure 5: Schematic of finite-volume computational cell in cylindrical geometry. Surface to volume ratios 
are independent of the wedge angle. The momentum source term in the axi-symmetric formulation of the Euler 
equations is proportional to dSt, and behaves as 1/Ar near the axis. 



Figure 7: Effect of viscosity for low nozzle pressures. The top section shows the results for a cutoff of 
gradients at a scale comparable to 0.1 mfp. The bottom section uses a cutoff at 10 mfp. The case on the left 
is for a nozzle pressure of 100 /iTorr, the case on the right is for 100 milliTorr (the experimental condition for 
the 16” shock tunnel). The cutoff allows greater stability and greater time steps in computing implicitely the 
viscous terms. The contours shown on the right hand side show that for our case of interest, the results are 
insensitive to the choice of cutoff. 
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Figure 6: Driven/Nozzle flow transient. This computation is obtained after coarsening the grid in the 
radial direction near the axis. The shock convergence on the axis (70 fisec) now leads to a stable shock structure 
(’bubble’) near the axis. The formation of an intense vortex is evident for 200 jzsec. Plots at later times would 
show that the reflected shock becomes essentially planar at ~ 500/isec, while the vortex structure persists. 




Figure 8: Nozzle flow transient (Mach contours). The curved primary shock can be identified on the right. 
Weak oblique shocks are emanating from the walls and converge to form a strong (upstream facing) mach disk 
on the axis, and reflect back. Strongly sheared flow is visible between the primary shock and the first oblique 
shocks. A contact discontinuity behind the primary shock has been numerically diffused beyond identification. 
A small vertical break in contours is due to an interpolation error by the graphics, between two different grids. 



Figure 9 Experimental shadowgraph, taken from [6], The upstream facing shock and the attached oblique 
shocks are clearly defined in this picture. Contact discontinuities are also very clear, and show strong vortices. 
All features of Figure 8 are contained in this picture. 
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t = 1.275 ms 
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Figure 10: Interaction in test chamber (Mach contours). In this example, as the primary shock diffracts at 
the end of the cone, the remainder of the shock structure has just reached the apex of the cone. For this case, 
steady flow over the testing body is achieved within 400 /^seconds. 
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Figure 11: Sensitivity of various flow quantities to effective nozzle area. The nozzle shape is replaced by 
a modified nozzle contour from the throat on, by removing a fixed fraction of the radius (i.e. the 87% nozzle 
removes 13% of the radius), to account for the boundary layer. Examination of the plots shows that the static 
pressure is the most sensitive variable, and therefore the most useful variable to determine experimentally. 
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ERATURE STATIC PRESSURE 
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Figure 12:Comparison of quasi- ID solutions with 2D axi-symmetric solution. Except for the flow within 
the boundary layer, the agreement is very good for the 87% case. Most variables have a flat profile in the 
core region, except the static pressure. Again, measurement of the static pressure should yield some valuable 
information on the nozzle flow, and for code validation. 





TEMPERATURE AT STATION N3 MOLE FRACTION OF NO AT STATION N3 



Figure 13. Temperature and N O mole fraction at two stations. There is still some significant core variation 
of the temperature at N3. The NO concentration profile remains essentially frozen, with most of the gradients 
in the boundary layer, as expected. 
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Figure 14: Computed Intensity at station N3, with and without the boundary layer. Figure 14-a (top), 
shows the intensity in absolute value, with the boundary layer included. Figure 14-b (bottom), shows the ratio 
of computed intensity with the boundary layer, versus the computed intensity without boundary layer. A similar 
plot for the exit plane would show a ratio at peak absorption of 0.2. 
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Figure 15: Comparison of quasi-ID and axi-symmetric cases (without BL). Figure 15-a (top) shows the ra- 
tio of the computed intensity for the quasi-ID solution (87%) versus the computed intensity for the 2D solution 
at station N3. Figure 15-b (bottom) shows a similar quantity for the exit plane. 
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Figure 17. Shock and Contact Discontinuity for same time sequence as Figure 16. 




Figure 18 : Comparison between experimental and computed pressure traces at two locations. Figure 18-a 
(top) shows the comparison at a few cm upstream of the main diaphragm in the driver tube. Figure 18-b (bot- 
tom) shows the traces for a point near the end of the driven tube. Significantly better agreement is obtained 
after assuming an ad-hoc pressure gradient within the driver tube, at the end of combustion. 
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